diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6fc1a8305..2d5913c06 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -101,20 +101,33 @@ message(STATUS "STORM_CPP_TESTS_BASE_PATH is ${STORM_CPP_TESTS_BASE_PATH}")
 
 # Main Sources
 file(GLOB_RECURSE STORM_HEADERS ${PROJECT_SOURCE_DIR}/src/*.h)
-file(GLOB_RECURSE STORM_SOURCES ${PROJECT_SOURCE_DIR}/src/*.cpp)
+file(GLOB_RECURSE STORM_SOURCES_WITHOUT_MAIN ${PROJECT_SOURCE_DIR}/src/*/*.cpp)
+file(GLOB_RECURSE STORM_MAIN_FILE ${PROJECT_SOURCE_DIR}/src/storm.cpp)
+set(STORM_SOURCES "${STORM_SOURCES_WITHOUT_MAIN};${STORM_MAIN_FILE}")
+file(GLOB_RECURSE STORM_ADAPTERS_FILES ${PROJECT_SOURCE_DIR}/src/adapters/*.h ${PROJECT_SOURCE_DIR}/src/adapters/*.cpp)
+file(GLOB_RECURSE STORM_EXCEPTIONS_FILES ${PROJECT_SOURCE_DIR}/src/exceptions/*.h ${PROJECT_SOURCE_DIR}/src/exceptions/*.cpp)
+file(GLOB_RECURSE STORM_FORMULA_FILES ${PROJECT_SOURCE_DIR}/src/formula/*.h ${PROJECT_SOURCE_DIR}/src/formula/*.cpp)
+file(GLOB_RECURSE STORM_MODELCHECKER_FILES ${PROJECT_SOURCE_DIR}/src/modelchecker/*.h ${PROJECT_SOURCE_DIR}/src/modelchecker/*.cpp)
+file(GLOB_RECURSE STORM_MODELS_FILES ${PROJECT_SOURCE_DIR}/src/models/*.h ${PROJECT_SOURCE_DIR}/src/models/*.cpp)
+file(GLOB_RECURSE STORM_PARSER_FILES ${PROJECT_SOURCE_DIR}/src/parser/*.h ${PROJECT_SOURCE_DIR}/src/parser/*.cpp)
+file(GLOB_RECURSE STORM_STORAGE_FILES ${PROJECT_SOURCE_DIR}/src/storage/*.h ${PROJECT_SOURCE_DIR}/src/storage/*.cpp)
+file(GLOB_RECURSE STORM_UTILITY_FILES ${PROJECT_SOURCE_DIR}/src/utility/*.h ${PROJECT_SOURCE_DIR}/src/utility/*.cpp)
 
 # Test Sources
 # Note that the tests also need the source files, except for the main file
-file(GLOB_RECURSE STORM_TEST_HEADERS ${PROJECT_SOURCE_DIR}/test/*.h)
-file(GLOB_RECURSE STORM_TEST_SOURCES ${PROJECT_SOURCE_DIR}/test/*.cpp ${PROJECT_SOURCE_DIR}/src/*/*.cpp)
-
-# Main Grouping
-source_group(Headers FILES ${STORM_HEADERS})
-source_group(Sources FILES ${STORM_SOURCES})
-
-# Test Grouping
-source_group(Headers FILES ${STORM_TEST_HEADERS})
-source_group(Sources FILES ${STORM_TEST_SOURCES})
+file(GLOB_RECURSE STORM_TEST_FILES ${PROJECT_SOURCE_DIR}/test/*.h ${PROJECT_SOURCE_DIR}/test/*.cpp)
+
+# Group the headers and sources
+source_group(main FILES ${STORM_MAIN_FILE})
+source_group(adapters FILES ${STORM_ADAPTERS_FILES})
+source_group(exceptions FILES ${STORM_EXCEPTIONS_FILES})
+source_group(formula FILES ${STORM_FORMULA_FILES})
+source_group(modelchecker FILES ${STORM_MODELCHECKER_FILES})
+source_group(models FILES ${STORM_MODELS_FILES})
+source_group(parser FILES ${STORM_PARSER_FILES})
+source_group(storage FILES ${STORM_STORAGE_FILES})
+source_group(utility FILES ${STORM_UTILITY_FILES})
+source_group(test FILES ${STORM_TEST_FILES})
 
 # Add base folder for better inclusion paths
 include_directories("${PROJECT_SOURCE_DIR}")
@@ -149,7 +162,7 @@ include_directories(${GMMXX_INCLUDE_DIR})
 # Add the executables
 # Must be created *after* Boost was added because of LINK_DIRECTORIES
 add_executable(storm ${STORM_SOURCES} ${STORM_HEADERS})
-add_executable(storm-tests ${STORM_TEST_SOURCES} ${STORM_TEST_HEADERS})
+add_executable(storm-tests ${STORM_TEST_FILES} ${STORM_SOURCES_WITHOUT_MAIN} ${STORM_HEADERS})
 
 # Add target link deps for Boost program options
 target_link_libraries(storm ${Boost_LIBRARIES})   
@@ -158,6 +171,10 @@ target_link_libraries(storm-tests ${Boost_LIBRARIES})
 set (STORM_USE_COTIRE ON)
 if (APPLE)
 	set(STORM_USE_COTIRE OFF)
+	
+    # Set up some Xcode specific settings
+	set(CMAKE_XCODE_ATTRIBUTE_CLANG_CXX_LANGUAGE_STANDARD "c++11")
+	set(CMAKE_XCODE_ATTRIBUTE_CLANG_CXX_LIBRARY "libc++")
 endif(APPLE)
 
 # Print Cotire Usage Status
diff --git a/src/formula/AbstractPathFormula.h b/src/formula/AbstractPathFormula.h
index d0ee74efd..459cbd057 100644
--- a/src/formula/AbstractPathFormula.h
+++ b/src/formula/AbstractPathFormula.h
@@ -61,7 +61,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T>* check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const = 0;
+	virtual std::vector<T>* check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const = 0;
 };
 
 } //namespace formula
diff --git a/src/formula/AbstractStateFormula.h b/src/formula/AbstractStateFormula.h
index 9947c5df6..cf15fdb30 100644
--- a/src/formula/AbstractStateFormula.h
+++ b/src/formula/AbstractStateFormula.h
@@ -58,7 +58,7 @@ public:
 	 *
 	 * @returns A bit vector indicating all states that satisfy the formula represented by the called object.
 	 */
-	virtual storm::storage::BitVector *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const = 0; // {
+	virtual storm::storage::BitVector *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const = 0; // {
 };
 
 } //namespace formula
diff --git a/src/formula/And.h b/src/formula/And.h
index f5d758c2b..80f46a821 100644
--- a/src/formula/And.h
+++ b/src/formula/And.h
@@ -162,7 +162,7 @@ public:
 	 *
 	 * @returns A bit vector indicating all states that satisfy the formula represented by the called object.
 	 */
-	virtual storm::storage::BitVector* check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const {
+	virtual storm::storage::BitVector* check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const {
 		return modelChecker.template as<IAndModelChecker>()->checkAnd(*this);
 	}
 	
diff --git a/src/formula/Ap.h b/src/formula/Ap.h
index 4950b4b72..0cb0d2480 100644
--- a/src/formula/Ap.h
+++ b/src/formula/Ap.h
@@ -98,7 +98,7 @@ public:
 	 *
 	 * @returns A bit vector indicating all states that satisfy the formula represented by the called object.
 	 */
-	virtual storm::storage::BitVector *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const {
+	virtual storm::storage::BitVector *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const {
 		return modelChecker.template as<IApModelChecker>()->checkAp(*this);
 	}
 	
diff --git a/src/formula/BoundedEventually.h b/src/formula/BoundedEventually.h
index fe7776f33..33c42f1ba 100644
--- a/src/formula/BoundedEventually.h
+++ b/src/formula/BoundedEventually.h
@@ -157,7 +157,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T>* check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T>* check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<IBoundedEventuallyModelChecker>()->checkBoundedEventually(*this, qualitative);
 	}
 	
diff --git a/src/formula/BoundedNaryUntil.h b/src/formula/BoundedNaryUntil.h
index 7c98af61e..9f83df869 100644
--- a/src/formula/BoundedNaryUntil.h
+++ b/src/formula/BoundedNaryUntil.h
@@ -179,7 +179,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<IBoundedNaryUntilModelChecker>()->checkBoundedNaryUntil(*this, qualitative);
 	}
 	
diff --git a/src/formula/BoundedUntil.h b/src/formula/BoundedUntil.h
index a5eec8576..685b1bfe1 100644
--- a/src/formula/BoundedUntil.h
+++ b/src/formula/BoundedUntil.h
@@ -184,7 +184,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<IBoundedUntilModelChecker>()->checkBoundedUntil(*this, qualitative);
 	}
 	
diff --git a/src/formula/CumulativeReward.h b/src/formula/CumulativeReward.h
index 3c38991a8..5d3dfa667 100644
--- a/src/formula/CumulativeReward.h
+++ b/src/formula/CumulativeReward.h
@@ -120,7 +120,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<ICumulativeRewardModelChecker>()->checkCumulativeReward(*this, qualitative);
 	}
 	
diff --git a/src/formula/Eventually.h b/src/formula/Eventually.h
index 668546f0c..8db4f1b01 100644
--- a/src/formula/Eventually.h
+++ b/src/formula/Eventually.h
@@ -131,7 +131,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<IEventuallyModelChecker>()->checkEventually(*this, qualitative);
 	}
 	
diff --git a/src/formula/Formulas.h b/src/formula/Formulas.h
index e5864bbc0..4ad64b1e7 100644
--- a/src/formula/Formulas.h
+++ b/src/formula/Formulas.h
@@ -38,7 +38,4 @@
 #include "AbstractPathFormula.h"
 #include "AbstractStateFormula.h"
 
-// FIXME: Why is this needed? To me this makes no sense.
-#include "modelchecker/DtmcPrctlModelChecker.h"
-
 #endif /* STORM_FORMULA_FORMULAS_H_ */
diff --git a/src/formula/Globally.h b/src/formula/Globally.h
index 7b58580e9..31d06bd05 100644
--- a/src/formula/Globally.h
+++ b/src/formula/Globally.h
@@ -132,7 +132,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<IGloballyModelChecker>()->checkGlobally(*this, qualitative);
 	}
 	
diff --git a/src/formula/InstantaneousReward.h b/src/formula/InstantaneousReward.h
index ed0f7affa..c2001540f 100644
--- a/src/formula/InstantaneousReward.h
+++ b/src/formula/InstantaneousReward.h
@@ -121,7 +121,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<IInstantaneousRewardModelChecker>()->checkInstantaneousReward(*this, qualitative);
 	}
 	
diff --git a/src/formula/Next.h b/src/formula/Next.h
index 91d68d05a..5f36b053e 100644
--- a/src/formula/Next.h
+++ b/src/formula/Next.h
@@ -133,7 +133,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<INextModelChecker>()->checkNext(*this, qualitative);
 	}
 	
diff --git a/src/formula/NoBoundOperator.h b/src/formula/NoBoundOperator.h
index 5261f3736..cbbc9f05f 100644
--- a/src/formula/NoBoundOperator.h
+++ b/src/formula/NoBoundOperator.h
@@ -136,7 +136,7 @@ public:
 	 *
 	 * @returns A vector indicating all states that satisfy the formula represented by the called object.
 	 */
-	virtual std::vector<T>* check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const {
+	virtual std::vector<T>* check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const {
 		return modelChecker.template as<INoBoundOperatorModelChecker>()->checkNoBoundOperator(*this);
 	}
 
diff --git a/src/formula/Not.h b/src/formula/Not.h
index 0ddb803c9..dde509ec8 100644
--- a/src/formula/Not.h
+++ b/src/formula/Not.h
@@ -127,7 +127,7 @@ public:
 	 *
 	 * @returns A bit vector indicating all states that satisfy the formula represented by the called object.
 	 */
-	virtual storm::storage::BitVector *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const {
+	virtual storm::storage::BitVector *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const {
 		return modelChecker.template as<INotModelChecker>()->checkNot(*this);  
 	}
 	
diff --git a/src/formula/Or.h b/src/formula/Or.h
index 52cc1d8a8..bb91bfa59 100644
--- a/src/formula/Or.h
+++ b/src/formula/Or.h
@@ -160,7 +160,7 @@ public:
 	 *
 	 * @returns A bit vector indicating all states that satisfy the formula represented by the called object.
 	 */
-	virtual storm::storage::BitVector *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const {
+	virtual storm::storage::BitVector *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const {
 		return modelChecker.template as<IOrModelChecker>()->checkOr(*this);
 	}
 	
diff --git a/src/formula/ProbabilisticBoundOperator.h b/src/formula/ProbabilisticBoundOperator.h
index 14429c677..8d651ef2c 100644
--- a/src/formula/ProbabilisticBoundOperator.h
+++ b/src/formula/ProbabilisticBoundOperator.h
@@ -117,7 +117,7 @@ public:
 	 *
 	 * @returns A bit vector indicating all states that satisfy the formula represented by the called object.
 	 */
-	virtual storm::storage::BitVector* check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const {
+	virtual storm::storage::BitVector* check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const {
 		return modelChecker.template as<IProbabilisticBoundOperatorModelChecker>()->checkProbabilisticBoundOperator(*this);
 	}
 };
diff --git a/src/formula/ReachabilityReward.h b/src/formula/ReachabilityReward.h
index ff8155eb5..f88deb30e 100644
--- a/src/formula/ReachabilityReward.h
+++ b/src/formula/ReachabilityReward.h
@@ -127,7 +127,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<IReachabilityRewardModelChecker>()->checkReachabilityReward(*this, qualitative);
 	}
 	
diff --git a/src/formula/RewardBoundOperator.h b/src/formula/RewardBoundOperator.h
index f602f89b9..18d715fbc 100644
--- a/src/formula/RewardBoundOperator.h
+++ b/src/formula/RewardBoundOperator.h
@@ -113,7 +113,7 @@ public:
 	 *
 	 * @returns A bit vector indicating all states that satisfy the formula represented by the called object.
 	 */
-	virtual storm::storage::BitVector* check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const {
+	virtual storm::storage::BitVector* check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const {
 		return modelChecker.template as<IRewardBoundOperatorModelChecker>()->checkRewardBoundOperator(*this);
 	}
 };
diff --git a/src/formula/StateBoundOperator.h b/src/formula/StateBoundOperator.h
index 9cfcfb862..323994c76 100644
--- a/src/formula/StateBoundOperator.h
+++ b/src/formula/StateBoundOperator.h
@@ -171,7 +171,7 @@ public:
 	 *
 	 * @returns A bit vector indicating all states that satisfy the formula represented by the called object.
 	 */
-	virtual storm::storage::BitVector *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const {
+	virtual storm::storage::BitVector *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const {
 		return modelChecker.template as<IStateBoundOperatorModelChecker>()->checkStateBoundOperator(*this);
 	}
 	
diff --git a/src/formula/SteadyStateOperator.h b/src/formula/SteadyStateOperator.h
index 755f12249..402a33c01 100644
--- a/src/formula/SteadyStateOperator.h
+++ b/src/formula/SteadyStateOperator.h
@@ -107,7 +107,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual storm::storage::BitVector* check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker) const {
+	virtual storm::storage::BitVector* check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker) const {
 		return modelChecker.template as<ISteadyStateOperatorModelChecker>()->checkSteadyStateOperator(*this);
 	}
 	
diff --git a/src/formula/SteadyStateReward.h b/src/formula/SteadyStateReward.h
index 799bc50a7..eae0def8b 100644
--- a/src/formula/SteadyStateReward.h
+++ b/src/formula/SteadyStateReward.h
@@ -84,7 +84,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<ISteadyStateRewardModelChecker>()->checkSteadyStateReward(*this, qualitative);
 	}
 
diff --git a/src/formula/Until.h b/src/formula/Until.h
index 562d397fc..b2efe6f45 100644
--- a/src/formula/Until.h
+++ b/src/formula/Until.h
@@ -158,7 +158,7 @@ public:
 	 *
 	 * @returns A vector indicating the probability that the formula holds for each state.
 	 */
-	virtual std::vector<T> *check(const storm::modelChecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
+	virtual std::vector<T> *check(const storm::modelchecker::AbstractModelChecker<T>& modelChecker, bool qualitative) const {
 		return modelChecker.template as<IUntilModelChecker>()->checkUntil(*this, qualitative);
 	}
 	
diff --git a/src/modelchecker/AbstractModelChecker.h b/src/modelchecker/AbstractModelChecker.h
index a2a41f1e3..de32619f7 100644
--- a/src/modelchecker/AbstractModelChecker.h
+++ b/src/modelchecker/AbstractModelChecker.h
@@ -1,5 +1,5 @@
 /*
- * DtmcPrctlModelChecker.h
+ * AbstractModelChecker.h
  *
  *  Created on: 22.10.2012
  *      Author: Thomas Heinemann
@@ -8,31 +8,40 @@
 #ifndef STORM_MODELCHECKER_ABSTRACTMODELCHECKER_H_
 #define STORM_MODELCHECKER_ABSTRACTMODELCHECKER_H_
 
-namespace storm { namespace modelChecker {
-template <class Type> class AbstractModelChecker;
-}}
+// Forward declaration of abstract model checker class needed by the formula classes.
+namespace storm {
+namespace modelchecker {
+	template <class Type> class AbstractModelChecker;
+}
+}
 
 #include "src/exceptions/InvalidPropertyException.h"
 #include "src/formula/Formulas.h"
 #include "src/storage/BitVector.h"
 #include "src/models/AbstractModel.h"
 
+#include "log4cplus/logger.h"
+#include "log4cplus/loggingmacros.h"
+
 #include <iostream>
 
+extern log4cplus::Logger logger;
+
 namespace storm {
-namespace modelChecker {
+namespace modelchecker {
 
 /*!
  * @brief
- * Interface for model checker classes.
+ * (Abstract) interface for all model checker classes.
  *
- * This class provides basic functions that are the same for all subclasses, but mainly only declares
- * abstract methods that are to be implemented in concrete instances.
- *
- * @attention This class is abstract.
+ * This class provides basic functions that are common to all model checkers (i.e. subclasses). It mainly declares
+ * abstract methods that are implemented in the concrete subclasses, but also covers checking procedures that are common
+ * to all model checkers for state-based models.
  */
 template<class Type>
 class AbstractModelChecker :
+	// A list of interfaces the model checker supports. Typically, for each of the interfaces, a check method needs to
+	// be implemented that performs the corresponding check.
 	public virtual storm::formula::IApModelChecker<Type>,
 	public virtual storm::formula::IAndModelChecker<Type>,
 	public virtual storm::formula::IOrModelChecker<Type>,
@@ -51,36 +60,70 @@ class AbstractModelChecker :
 	public virtual storm::formula::IInstantaneousRewardModelChecker<Type> {
 	
 public:
-	explicit AbstractModelChecker(storm::models::AbstractModel<Type>& model)
-		: model(model) {
-		// Nothing to do here...
+	/*!
+	 * Constructs an AbstractModelChecker with the given model.
+	 */
+	explicit AbstractModelChecker(storm::models::AbstractModel<Type> const& model) : model(model) {
+		// Intentionally left empty.
 	}
 	
-	explicit AbstractModelChecker(AbstractModelChecker<Type>* modelChecker)
-		: model(modelChecker->model) {
+	/*!
+	 * Copy constructs an AbstractModelChecker from the given model checker. In particular, this means that the newly
+	 * constructed model checker will have the model of the given model checker as its associated model.
+	 */
+	explicit AbstractModelChecker(AbstractModelChecker<Type> const& modelchecker) : model(modelchecker.model) {
+		// Intentionally left empty.
 	}
 	
+	/*!
+	 * Virtual destructor. Needs to be virtual, because this class has virtual methods.
+	 */
 	virtual ~AbstractModelChecker() {
-		//intentionally left empty
+		// Intentionally left empty.
 	}
 
+	/*!
+	 * 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.
+	 * 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>
 	const Target<Type>* as() const {
 		try {
 			const Target<Type>* target = dynamic_cast<const Target<Type>*>(this);
 			return target;
 		} catch (std::bad_cast& bc) {
-			std::cerr << "Bad cast: tried to cast " << typeid(*this).name() << " to " << typeid(Target<Type>).name() << std::endl;
+			LOG4CPLUS_ERROR(logger, "Bad cast: tried to cast " << typeid(*this).name() << " to " << typeid(Target<Type>).name() << ".");
+			throw bc;
 		}
 		return nullptr;
 	}
 
 	/*!
-	 * Checks the given state formula on the DTMC and prints the result (true/false) for all initial
-	 * states.
+	 * 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
+	 * is not of the requested type, type casting will fail and result in an exception.
+	 */
+	template <class Model>
+	Model const& getModel() const {
+		try {
+			Model const& target = dynamic_cast<Model const&>(this->model);
+			return target;
+		} catch (std::bad_cast& bc) {
+			LOG4CPLUS_ERROR(logger, "Bad cast: tried to cast " << typeid(this->model).name() << " to " << typeid(Model).name() << ".");
+			throw bc;
+		}
+	}
+
+	/*!
+	 * Checks the given state formula on the model and prints the result (true/false) for all initial states, i.e.
+	 * states that carry the atomic proposition "init".
+	 *
 	 * @param stateFormula The formula to be checked.
 	 */
-	void check(const storm::formula::AbstractStateFormula<Type>& stateFormula) const {
+	void check(storm::formula::AbstractStateFormula<Type> const& stateFormula) const {
 		std::cout << std::endl;
 		LOG4CPLUS_INFO(logger, "Model checking formula\t" << stateFormula.toString());
 		std::cout << "Model checking formula:\t" << stateFormula.toString() << std::endl;
@@ -96,6 +139,7 @@ public:
 			delete result;
 		} catch (std::exception& e) {
 			std::cout << "Error during computation: " << e.what() << "Skipping property." << std::endl;
+			LOG4CPLUS_ERROR(logger, "Error during computation: " << e.what() << "Skipping property.");
 			if (result != nullptr) {
 				delete result;
 			}
@@ -105,11 +149,12 @@ public:
 	}
 
 	/*!
-	 * Checks the given operator (with no bound) on the DTMC and prints the result
-	 * (probability/rewards) for all initial states.
+	 * Checks the given formula (with no bound) on the model and prints the result (probability/rewards) for all
+	 * initial states, i.e. states that carry the atomic proposition "init".
+	 *
 	 * @param noBoundFormula The formula to be checked.
 	 */
-	void check(const storm::formula::NoBoundOperator<Type>& noBoundFormula) const {
+	void check(storm::formula::NoBoundOperator<Type> const& noBoundFormula) const {
 		std::cout << std::endl;
 		LOG4CPLUS_INFO(logger, "Model checking formula\t" << noBoundFormula.toString());
 		std::cout << "Model checking formula:\t" << noBoundFormula.toString() << std::endl;
@@ -134,15 +179,15 @@ public:
 	}
 
 	/*!
-	 * The check method for a formula with an AP node as root in its formula tree
+	 * Checks the given formula consisting of a single atomic proposition.
 	 *
-	 * @param formula The Ap state formula to check
-	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 * @param formula The formula to be checked.
+	 * @returns The set of states satisfying the formula represented by a bit vector.
 	 */
-	storm::storage::BitVector* checkAp(const storm::formula::Ap<Type>& formula) const {
-		if (formula.getAp().compare("true") == 0) {
+	storm::storage::BitVector* checkAp(storm::formula::Ap<Type> const& formula) const {
+		if (formula.getAp() == "true") {
 			return new storm::storage::BitVector(model.getNumberOfStates(), true);
-		} else if (formula.getAp().compare("false") == 0) {
+		} else if (formula.getAp() == "false") {
 			return new storm::storage::BitVector(model.getNumberOfStates());
 		}
 
@@ -155,12 +200,12 @@ public:
 	}
 
 	/*!
-	 * The check method for a state formula with an And node as root in its formula tree
+	 * Checks the given formula that is a logical "and" of two formulae.
 	 *
-	 * @param formula The And formula to check
-	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 * @param formula The formula to be checked.
+	 * @returns The set of states satisfying the formula represented by a bit vector.
 	 */
-	storm::storage::BitVector* checkAnd(const storm::formula::And<Type>& formula) const {
+	storm::storage::BitVector* checkAnd(storm::formula::And<Type> const& formula) const {
 		storm::storage::BitVector* result = formula.getLeft().check(*this);
 		storm::storage::BitVector* right = formula.getRight().check(*this);
 		(*result) &= (*right);
@@ -169,43 +214,43 @@ public:
 	}
 
 	/*!
-	 * The check method for a formula with a Not node as root in its formula tree
+	 * Checks the given formula that is a logical "or" of two formulae.
 	 *
-	 * @param formula The Not state formula to check
-	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 * @param formula The formula to check.
+	 * @returns The set of states satisfying the formula represented by a bit vector.
 	 */
-	storm::storage::BitVector* checkNot(const storm::formula::Not<Type>& formula) const {
-		storm::storage::BitVector* result = formula.getChild().check(*this);
-		result->complement();
+	virtual storm::storage::BitVector* checkOr(storm::formula::Or<Type> const& formula) const {
+		storm::storage::BitVector* result = formula.getLeft().check(*this);
+		storm::storage::BitVector* right = formula.getRight().check(*this);
+		(*result) |= (*right);
+		delete right;
 		return result;
 	}
 
 	/*!
-	 * The check method for a state formula with an Or node as root in its formula tree
+	 * Checks the given formula that is a logical "not" of a sub-formula.
 	 *
-	 * @param formula The Or state formula to check
-	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 * @param formula The formula to check.
+	 * @returns The set of states satisfying the formula represented by a bit vector.
 	 */
-	virtual storm::storage::BitVector* checkOr(const storm::formula::Or<Type>& formula) const {
-		storm::storage::BitVector* result = formula.getLeft().check(*this);
-		storm::storage::BitVector* right = formula.getRight().check(*this);
-		(*result) |= (*right);
-		delete right;
+	storm::storage::BitVector* checkNot(const storm::formula::Not<Type>& formula) const {
+		storm::storage::BitVector* result = formula.getChild().check(*this);
+		result->complement();
 		return result;
 	}
 
+
 	/*!
-	 * The check method for a state formula with a bound operator node as root in
-	 * its formula tree
+	 * Checks the given formula that is a P operator over a path formula featuring a value bound.
 	 *
-	 * @param formula The state formula to check
-	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 * @param formula The formula to check.
+	 * @returns The set of states satisfying the formula represented by a bit vector.
 	 */
-	storm::storage::BitVector* checkProbabilisticBoundOperator(const storm::formula::ProbabilisticBoundOperator<Type>& formula) const {
+	storm::storage::BitVector* checkProbabilisticBoundOperator(storm::formula::ProbabilisticBoundOperator<Type> const& formula) const {
 		// First, we need to compute the probability for satisfying the path formula for each state.
 		std::vector<Type>* quantitativeResult = formula.getPathFormula().check(*this, false);
 
-		// Create resulting bit vector, which will hold the yes/no-answer for every state.
+		// Create resulting bit vector that will hold the yes/no-answer for every state.
 		storm::storage::BitVector* result = new storm::storage::BitVector(quantitativeResult->size());
 
 		// Now, we can compute which states meet the bound specified in this operator and set the
@@ -222,17 +267,16 @@ public:
 	}
 
 	/*!
-	 * The check method for a state formula with a bound operator node as root in
-	 * its formula tree
+	 * Checks the given formula that is an R operator over a reward formula featuring a value bound.
 	 *
-	 * @param formula The state formula to check
-	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 * @param formula The formula to check.
+	 * @returns The set of states satisfying the formula represented by a bit vector.
 	 */
 	storm::storage::BitVector* checkRewardBoundOperator(const storm::formula::RewardBoundOperator<Type>& formula) const {
 		// First, we need to compute the probability for satisfying the path formula for each state.
 		std::vector<Type>* quantitativeResult = formula.getPathFormula().check(*this, false);
 
-		// Create resulting bit vector, which will hold the yes/no-answer for every state.
+		// Create resulting bit vector that will hold the yes/no-answer for every state.
 		storm::storage::BitVector* result = new storm::storage::BitVector(quantitativeResult->size());
 
 		// Now, we can compute which states meet the bound specified in this operator and set the
@@ -248,22 +292,18 @@ public:
 		return result;
 	}
 	
-	void setModel(storm::models::AbstractModel<Type>& model) {
-		this->model = model;
-	}
-	
-	template <class Model>
-	Model& getModel() const {
-		return *dynamic_cast<Model*>(&this->model);
-	}
-
 private:
-	storm::models::AbstractModel<Type>& model;
 
+	/*!
+	 * A constant reference to the model associated with this model checker.
+	 *
+	 * @note that we do not own this object, but merely have a constant reference to  it. That means that using the
+	 * model checker object is unsafe after the object has been destroyed.
+	 */
+	storm::models::AbstractModel<Type> const& model;
 };
 
-} //namespace modelChecker
-
-} //namespace storm
+} // namespace modelchecker
+} // namespace storm
 
 #endif /* STORM_MODELCHECKER_DTMCPRCTLMODELCHECKER_H_ */
diff --git a/src/modelchecker/DtmcPrctlModelChecker.h b/src/modelchecker/DtmcPrctlModelChecker.h
deleted file mode 100644
index 3fa4cc238..000000000
--- a/src/modelchecker/DtmcPrctlModelChecker.h
+++ /dev/null
@@ -1,411 +0,0 @@
-/*
- * DtmcPrctlModelChecker.h
- *
- *  Created on: 22.10.2012
- *      Author: Thomas Heinemann
- */
-
-#ifndef STORM_MODELCHECKER_DTMCPRCTLMODELCHECKER_H_
-#define STORM_MODELCHECKER_DTMCPRCTLMODELCHECKER_H_
-
-#include <vector>
-
-#include "src/formula/Formulas.h"
-#include "src/utility/Vector.h"
-#include "src/storage/SparseMatrix.h"
-
-#include "src/models/Dtmc.h"
-#include "src/storage/BitVector.h"
-#include "src/exceptions/InvalidPropertyException.h"
-#include "src/utility/Vector.h"
-#include "src/utility/GraphAnalyzer.h"
-#include "src/modelchecker/AbstractModelChecker.h"
-
-#include "log4cplus/logger.h"
-#include "log4cplus/loggingmacros.h"
-extern log4cplus::Logger logger;
-
-namespace storm {
-
-namespace modelChecker {
-
-/*!
- * @brief
- * Interface for model checker classes.
- *
- * This class provides basic functions that are the same for all subclasses, but mainly only declares
- * abstract methods that are to be implemented in concrete instances.
- *
- * @attention This class is abstract.
- */
-template<class Type>
-class DtmcPrctlModelChecker : public AbstractModelChecker<Type> {
-
-public:
-	/*!
-	 * Constructor
-	 *
-	 * @param model The dtmc model which is checked.
-	 */
-	explicit DtmcPrctlModelChecker(storm::models::Dtmc<Type>& model) : AbstractModelChecker<Type>(model) {
-		// Intentionally left empty.
-	}
-
-	/*!
-	 * Copy constructor
-	 *
-	 * @param modelChecker The model checker that is copied.
-	 */
-	explicit DtmcPrctlModelChecker(const storm::modelChecker::DtmcPrctlModelChecker<Type>* modelChecker) : AbstractModelChecker<Type>(modelChecker) {
-		// Intentionally left empty.
-	}
-
-	/*!
-	 * Destructor
-	 */
-	virtual ~DtmcPrctlModelChecker() {
-		// Intentionally left empty.
-	}
-
-	/*!
-	 * @returns A reference to the dtmc of the model checker.
-	 */
-	storm::models::Dtmc<Type>& getModel() const {
-		return AbstractModelChecker<Type>::template getModel<storm::models::Dtmc<Type>>();
-	}
-
-	/*!
-	 * The check method for a state formula with a probabilistic operator node without bounds as root
-	 * in its formula tree
-	 *
-	 * @param formula The state formula to check
-	 * @returns The set of states satisfying the formula, represented by a bit vector
-	 */
-	std::vector<Type>* checkNoBoundOperator(const storm::formula::NoBoundOperator<Type>& formula) const {
-		// Check if the operator was an optimality operator and report a warning in that case.
-		if (formula.isOptimalityOperator()) {
-			LOG4CPLUS_WARN(logger, "Formula contains min/max operator which is not meaningful over deterministic models.");
-		}
-		return formula.getPathFormula().check(*this, false);
-	}
-
-	/*!
-	 * The check method for a path formula with a Bounded Until operator node as root in its formula tree
-	 *
-	 * @param formula The Bounded Until path formula to check
-	 * @returns for each state the probability that the path formula holds.
-	 */
-	virtual std::vector<Type>* checkBoundedUntil(const storm::formula::BoundedUntil<Type>& formula, bool qualitative) const {
-		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
-		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
-		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
-
-		// Copy the matrix before we make any changes.
-		storm::storage::SparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionMatrix());
-
-		// Make all rows absorbing that violate both sub-formulas or satisfy the second sub-formula.
-		tmpMatrix.makeRowsAbsorbing(~(*leftStates | *rightStates) | *rightStates);
-
-		// Delete obsolete intermediates.
-		delete leftStates;
-		delete rightStates;
-
-		// Create the vector with which to multiply.
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
-		storm::utility::setVectorValues(result, *rightStates, storm::utility::constGetOne<Type>());
-
-		// Perform the matrix vector multiplication as often as required by the formula bound.
-		this->performMatrixVectorMultiplication(tmpMatrix, *result, nullptr, formula.getBound());
-
-		// Return result.
-		return result;
-	}
-
-	/*!
-	 * The check method for a path formula with a Next operator node as root in its formula tree
-	 *
-	 * @param formula The Next path formula to check
-	 * @returns for each state the probability that the path formula holds.
-	 */
-	virtual std::vector<Type>* checkNext(const storm::formula::Next<Type>& formula, bool qualitative) const {
-		// First, we need to compute the states that satisfy the child formula of the next-formula.
-		storm::storage::BitVector* nextStates = formula.getChild().check(*this);
-
-		// Create the vector with which to multiply and initialize it correctly.
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
-		storm::utility::setVectorValues(result, *nextStates, storm::utility::constGetOne<Type>());
-
-		// Delete obsolete intermediate.
-		delete nextStates;
-
-		// Perform one single matrix-vector multiplication.
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result);
-
-		// Return result.
-		return result;
-	}
-
-	/*!
-	 * The check method for a path formula with a Bounded Eventually operator node as root in its
-	 * formula tree
-	 *
-	 * @param formula The Bounded Eventually path formula to check
-	 * @returns for each state the probability that the path formula holds
-	 */
-	virtual std::vector<Type>* checkBoundedEventually(const storm::formula::BoundedEventually<Type>& formula, bool qualitative) const {
-		// Create equivalent temporary bounded until formula and check it.
-		storm::formula::BoundedUntil<Type> temporaryBoundedUntilFormula(new storm::formula::Ap<Type>("true"), formula.getChild().clone(), formula.getBound());
-		return this->checkBoundedUntil(temporaryBoundedUntilFormula, qualitative);
-	}
-
-	/*!
-	 * The check method for a path formula with an Eventually operator node as root in its formula tree
-	 *
-	 * @param formula The Eventually path formula to check
-	 * @returns for each state the probability that the path formula holds
-	 */
-	virtual std::vector<Type>* checkEventually(const storm::formula::Eventually<Type>& formula, bool qualitative) const {
-		// Create equivalent temporary until formula and check it.
-		storm::formula::Until<Type> temporaryUntilFormula(new storm::formula::Ap<Type>("true"), formula.getChild().clone());
-		return this->checkUntil(temporaryUntilFormula, qualitative);
-	}
-
-	/*!
-	 * The check method for a path formula with a Globally operator node as root in its formula tree
-	 *
-	 * @param formula The Globally path formula to check
-	 * @returns for each state the probability that the path formula holds
-	 */
-	virtual std::vector<Type>* checkGlobally(const storm::formula::Globally<Type>& formula, bool qualitative) const {
-		// Create "equivalent" temporary eventually formula and check it.
-		storm::formula::Eventually<Type> temporaryEventuallyFormula(new storm::formula::Not<Type>(formula.getChild().clone()));
-		std::vector<Type>* result = this->checkEventually(temporaryEventuallyFormula, qualitative);
-
-		// Now subtract the resulting vector from the constant one vector to obtain final result.
-		storm::utility::subtractFromConstantOneVector(result);
-		return result;
-	}
-
-	/*!
-	 * The check method for a path formula with an Until operator node as root in its formula tree
-	 *
-	 * @param formula The Until path formula to check
-	 * @returns for each state the probability that the path formula holds.
-	 */
-	virtual std::vector<Type>* checkUntil(const storm::formula::Until<Type>& formula, bool qualitative) const {
-		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
-		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
-		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
-
-		// Then, we need to identify the states which have to be taken out of the matrix, i.e.
-		// all states that have probability 0 and 1 of satisfying the until-formula.
-		storm::storage::BitVector statesWithProbability0(this->getModel().getNumberOfStates());
-		storm::storage::BitVector statesWithProbability1(this->getModel().getNumberOfStates());
-		storm::utility::GraphAnalyzer::performProb01(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
-
-		// Delete intermediate results that are obsolete now.
-		delete leftStates;
-		delete rightStates;
-
-		// Perform some logging.
-		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(maybeStatesSetBitCount);
-			this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, statesWithProbability1, &b);
-
-			this->solveEquationSystem(*submatrix, x, b);
-
-			// Delete the created submatrix.
-			delete submatrix;
-
-			// Set values of resulting vector according to result.
-			storm::utility::setVectorValues<Type>(result, maybeStates, x);
-		} 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::setVectorValues<Type>(result, maybeStates, Type(0.5));
-		}
-
-		// Set values of resulting vector that are known exactly.
-		storm::utility::setVectorValues<Type>(result, statesWithProbability0, storm::utility::constGetZero<Type>());
-		storm::utility::setVectorValues<Type>(result, statesWithProbability1, storm::utility::constGetOne<Type>());
-
-		return result;
-	}
-
-	/*!
-	 * The check method for a path formula with an Instantaneous Reward operator node as root in its
-	 * formula tree
-	 *
-	 * @param formula The Instantaneous Reward formula to check
-	 * @returns for each state the reward that the instantaneous reward yields
-	 */
-	virtual std::vector<Type>* checkInstantaneousReward(const storm::formula::InstantaneousReward<Type>& formula, bool qualitative) const {
-		// Only compute the result if the model has a state-based reward model.
-		if (!this->getModel().hasStateRewards()) {
-			LOG4CPLUS_ERROR(logger, "Missing (state-based) reward model for formula.");
-			throw storm::exceptions::InvalidPropertyException() << "Missing (state-based) reward model for formula.";
-		}
-
-		// Initialize result to state rewards of the model.
-		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
-
-		// 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());
-
-		// Return result.
-		return result;
-	}
-
-	/*!
-	 * The check method for a path formula with a Cumulative Reward operator node as root in its
-	 * formula tree
-	 *
-	 * @param formula The Cumulative Reward formula to check
-	 * @returns for each state the reward that the cumulative reward yields
-	 */
-	virtual std::vector<Type>* checkCumulativeReward(const storm::formula::CumulativeReward<Type>& formula, bool qualitative) const {
-		// Only compute the result if the model has at least one reward model.
-		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
-			LOG4CPLUS_ERROR(logger, "Missing reward model for formula.");
-			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
-		}
-
-		// Compute the reward vector to add in each step based on the available reward models.
-		std::vector<Type>* totalRewardVector = nullptr;
-		if (this->getModel().hasTransitionRewards()) {
-			totalRewardVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
-			if (this->getModel().hasStateRewards()) {
-				gmm::add(*this->getModel().getStateRewardVector(), *totalRewardVector);
-			}
-		} else {
-			totalRewardVector = new std::vector<Type>(*this->getModel().getStateRewardVector());
-		}
-
-		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
-
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, totalRewardVector, formula.getBound());
-
-		// Delete temporary variables and return result.
-		delete totalRewardVector;
-		return result;
-	}
-
-	/*!
-	 * The check method for a path formula with a Reachability Reward operator node as root in its
-	 * formula tree
-	 *
-	 * @param formula The Reachbility Reward formula to check
-	 * @returns for each state the reward that the reachability reward yields
-	 */
-	virtual std::vector<Type>* checkReachabilityReward(const storm::formula::ReachabilityReward<Type>& formula, bool qualitative) const {
-		// Only compute the result if the model has at least one reward model.
-		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
-			LOG4CPLUS_ERROR(logger, "Missing reward model for formula. Skipping formula");
-			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
-		}
-
-		// Determine the states for which the target predicate holds.
-		storm::storage::BitVector* targetStates = formula.getChild().check(*this);
-
-		// Determine which states have a reward of infinity by definition.
-		storm::storage::BitVector infinityStates(this->getModel().getNumberOfStates());
-		storm::storage::BitVector trueStates(this->getModel().getNumberOfStates(), true);
-		storm::utility::GraphAnalyzer::performProb1(this->getModel(), trueStates, *targetStates, &infinityStates);
-		infinityStates.complement();
-
-		// 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::selectVectorValues(&b, maybeStates, *pointwiseProductRowSumVector);
-				delete pointwiseProductRowSumVector;
-
-				if (this->getModel().hasStateRewards()) {
-					// If a state-based reward model is also available, we need to add this vector
-					// as well. As the state reward vector contains entries not just for the states
-					// that we still consider (i.e. maybeStates), we need to extract these values
-					// first.
-					std::vector<Type> subStateRewards(maybeStatesSetBitCount);
-					storm::utility::selectVectorValues(&subStateRewards, maybeStates, *this->getModel().getStateRewardVector());
-					gmm::add(subStateRewards, b);
-				}
-			} else {
-				// If only a state-based reward model is  available, we take this vector as the
-				// right-hand side. As the state reward vector contains entries not just for the
-				// states that we still consider (i.e. maybeStates), we need to extract these values
-				// first.
-				storm::utility::selectVectorValues(&b, maybeStates, *this->getModel().getStateRewardVector());
-			}
-
-			this->solveEquationSystem(*submatrix, x, b);
-
-			// Set values of resulting vector according to result.
-			storm::utility::setVectorValues<Type>(result, maybeStates, x);
-
-			// Delete temporary matrix and right-hand side.
-			delete submatrix;
-		}
-
-		// Set values of resulting vector that are known exactly.
-		storm::utility::setVectorValues(result, *targetStates, storm::utility::constGetZero<Type>());
-		storm::utility::setVectorValues(result, infinityStates, storm::utility::constGetInfinity<Type>());
-
-		// Delete temporary storages and return result.
-		delete targetStates;
-		return result;
-	}
-
-private:
-	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type>* summand = nullptr, uint_fast64_t repetitions = 1) const = 0;
-
-	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type> const& b) const = 0;
-};
-
-} //namespace modelChecker
-
-} //namespace storm
-
-#endif /* STORM_MODELCHECKER_DTMCPRCTLMODELCHECKER_H_ */
diff --git a/src/modelchecker/EigenDtmcPrctlModelChecker.h b/src/modelchecker/EigenDtmcPrctlModelChecker.h
index 167b02f72..311b69448 100644
--- a/src/modelchecker/EigenDtmcPrctlModelChecker.h
+++ b/src/modelchecker/EigenDtmcPrctlModelChecker.h
@@ -11,7 +11,7 @@
 #include "src/utility/Vector.h"
 
 #include "src/models/Dtmc.h"
-#include "src/modelchecker/DtmcPrctlModelChecker.h"
+#include "src/modelchecker/SparseDtmcPrctlModelChecker.h"
 #include "src/utility/GraphAnalyzer.h"
 #include "src/utility/ConstTemplates.h"
 #include "src/exceptions/NoConvergenceException.h"
@@ -29,20 +29,19 @@
 extern log4cplus::Logger logger;
 
 namespace storm {
-
-namespace modelChecker {
+namespace modelchecker {
 
 /*
  * A model checking engine that makes use of the eigen backend.
  */
 template <class Type>
-class EigenDtmcPrctlModelChecker : public DtmcPrctlModelChecker<Type> {
+class EigenDtmcPrctlModelChecker : public SparseDtmcPrctlModelChecker<Type> {
 
 typedef Eigen::Matrix<Type, -1, 1, 0, -1, 1> VectorType;
 typedef Eigen::Map<VectorType> MapType;
 
 public:
-	explicit EigenDtmcPrctlModelChecker(storm::models::Dtmc<Type>& dtmc) : DtmcPrctlModelChecker<Type>(dtmc) {
+	explicit EigenDtmcPrctlModelChecker(storm::models::Dtmc<Type>& dtmc) : SparseDtmcPrctlModelChecker<Type>(dtmc) {
 		// Intentionally left empty.
 	}
 
@@ -136,8 +135,7 @@ private:
 	}
 };
 
-} //namespace modelChecker
-
+} //namespace modelchecker
 } //namespace storm
 
 #endif /* STORM_MODELCHECKER_EIGENDTMCPRCTLMODELCHECKER_H_ */
diff --git a/src/modelchecker/ForwardDeclarations.h b/src/modelchecker/ForwardDeclarations.h
index 8239e2023..014d5eec0 100644
--- a/src/modelchecker/ForwardDeclarations.h
+++ b/src/modelchecker/ForwardDeclarations.h
@@ -2,25 +2,21 @@
  * ForwardDeclarations.h
  *
  *  Created on: 14.01.2013
- *      Author: thomas
+ *      Author: Thomas Heinemann
  */
 
 #ifndef STORM_MODELCHECKER_FORWARDDECLARATIONS_H_
 #define STORM_MODELCHECKER_FORWARDDECLARATIONS_H_
 
-
+// Forward declare the abstract model checker. This is used by the formula classes that need this declaration for
+// the callback methods (i.e., the check methods).
 namespace storm {
-
-namespace modelChecker {
+namespace modelchecker {
 
 template <class Type>
 class AbstractModelChecker;
 
-template <class Type>
-class DtmcPrctlModelChecker;
-
-} //namespace modelChecker
-
+} //namespace modelchecker
 } //namespace storm
 
 #endif /* STORM_MODELCHECKER_FORWARDDECLARATIONS_H_ */
diff --git a/src/modelchecker/GmmxxDtmcPrctlModelChecker.h b/src/modelchecker/GmmxxDtmcPrctlModelChecker.h
index 3c912a3c3..d07116929 100644
--- a/src/modelchecker/GmmxxDtmcPrctlModelChecker.h
+++ b/src/modelchecker/GmmxxDtmcPrctlModelChecker.h
@@ -8,48 +8,54 @@
 #ifndef STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_
 #define STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_
 
-#include <cmath>
-
-#include "src/models/Dtmc.h"
-#include "src/modelchecker/DtmcPrctlModelChecker.h"
-#include "src/utility/Vector.h"
-#include "src/utility/ConstTemplates.h"
-#include "src/utility/Settings.h"
+#include "src/modelchecker/SparseDtmcPrctlModelChecker.h"
 #include "src/adapters/GmmxxAdapter.h"
-#include "src/exceptions/InvalidPropertyException.h"
 #include "src/storage/JacobiDecomposition.h"
+#include "src/utility/ConstTemplates.h"
+#include "src/utility/Settings.h"
 
 #include "gmm/gmm_matrix.h"
 #include "gmm/gmm_iter_solvers.h"
 
-#include "log4cplus/logger.h"
-#include "log4cplus/loggingmacros.h"
-
-extern log4cplus::Logger logger;
+#include <cmath>
 
 namespace storm {
-
-namespace modelChecker {
+namespace modelchecker {
 
 /*
- * A model checking engine that makes use of the gmm++ backend.
+ * An implementation of the SparseDtmcPrctlModelChecker interface that uses gmm++ as the solving backend.
  */
 template <class Type>
-class GmmxxDtmcPrctlModelChecker
-	: public DtmcPrctlModelChecker<Type> {
+class GmmxxDtmcPrctlModelChecker : public SparseDtmcPrctlModelChecker<Type> {
 
 public:
-	explicit GmmxxDtmcPrctlModelChecker(storm::models::Dtmc<Type>& dtmc) : DtmcPrctlModelChecker<Type>(dtmc) {
+	/*!
+	 * Constructs a GmmxxDtmcPrctlModelChecker with the given model.
+	 *
+	 * @param model The DTMC to be checked.
+	 */
+	explicit GmmxxDtmcPrctlModelChecker(storm::models::Dtmc<Type>& dtmc) : SparseDtmcPrctlModelChecker<Type>(dtmc) {
 		// Intentionally left empty.
 	}
 
+	/*!
+	 * Copy constructs a GmmxxDtmcPrctlModelChecker from the given model checker. In particular, this means that the newly
+	 * constructed model checker will have the model of the given model checker as its associated model.
+	 */
+	explicit GmmxxDtmcPrctlModelChecker(storm::modelchecker::GmmxxDtmcPrctlModelChecker<Type> const& modelchecker) : SparseDtmcPrctlModelChecker<Type>(modelchecker) {
+		// Intentionally left empty.
+	}
+
+	/*!
+	 * Virtual destructor. Needs to be virtual, because this class has virtual methods.
+	 */
 	virtual ~GmmxxDtmcPrctlModelChecker() {
-		// Intentionally left empty
+		// Intentionally left empty.
 	}
 
 	/*!
 	 * Returns the name of this module.
-	 * @return The name of this module.
+	 * @returns The name of this module.
 	 */
 	static std::string getModuleName() {
 		return "gmm++det";
@@ -58,7 +64,7 @@ public:
 	/*!
 	 * Returns a trigger such that if the option "matrixlib" is set to "gmm++", this model checker
 	 * is to be used.
-	 * @return An option trigger for this module.
+	 * @returns An option trigger for this module.
 	 */
 	static std::pair<std::string, std::string> getOptionTrigger() {
 		return std::pair<std::string, std::string>("matrixlib", "gmm++");
@@ -97,28 +103,33 @@ public:
 
 private:
 
-	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type>* summand, uint_fast64_t repetitions = 1) const {
-		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(matrix);
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b, uint_fast64_t n = 1) const {
+		// Transform the transition probability A to the gmm++ format to use its arithmetic.
+		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
 
-		// Now perform matrix-vector multiplication as long as we meet the bound.
+		// Set up some temporary variables so that we can just swap pointers instead of copying the result after each
+		// iteration.
 		std::vector<Type>* swap = nullptr;
-		std::vector<Type>* currentVector = &vector;
+		std::vector<Type>* currentVector = &x;
 		std::vector<Type>* tmpVector = new std::vector<Type>(this->getModel().getNumberOfStates());
-		for (uint_fast64_t i = 0; i < repetitions; ++i) {
+
+		// Now perform matrix-vector multiplication as long as we meet the bound.
+		for (uint_fast64_t i = 0; i < n; ++i) {
 			gmm::mult(*gmmxxMatrix, *currentVector, *tmpVector);
 			swap = tmpVector;
 			tmpVector = currentVector;
 			currentVector = swap;
 
 			// If requested, add an offset to the current result vector.
-			if (summand != nullptr) {
-				gmm::add(*summand, *currentVector);
+			if (b != nullptr) {
+				gmm::add(*b, *currentVector);
 			}
 		}
 
-		if (repetitions % 2 == 1) {
-			std::swap(vector, *currentVector);
+		// If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, because
+		// the output is supposed to be stored in x.
+		if (n % 2 == 1) {
+			std::swap(x, *currentVector);
 			delete currentVector;
 		} else {
 			delete tmpVector;
@@ -127,17 +138,7 @@ private:
 		delete gmmxxMatrix;
 	}
 
-	/*!
-	 * Solves the linear equation system Ax=b with the given parameters.
-	 *
-	 * @param A The matrix A specifying the coefficients of the linear equations.
-	 * @param x The vector x for which to solve the equations. The initial value of the elements of
-	 * this vector are used as the initial guess and might thus influence performance and convergence.
-	 * @param b The vector b specifying the values on the right-hand-sides of the equations.
-	 * @return The solution of the system of linear equations in form of the elements of the vector
-	 * x.
-	 */
-	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type> const& b) const {
+	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const {
 		// Get the settings object to customize linear solving.
 		storm::settings::Settings* s = storm::settings::instance();
 
@@ -145,31 +146,38 @@ private:
 		// and the like.
 		gmm::iteration iter(s->get<double>("precision"), 0, s->get<unsigned>("maxiter"));
 
-		// Now do the actual solving.
-		LOG4CPLUS_INFO(logger, "Starting iterative solver.");
+		// Print some information about the used preconditioner.
 		const std::string& precond = s->getString("precond");
-		if (precond == "ilu") {
-			LOG4CPLUS_INFO(logger, "Using ILU preconditioner.");
-		} else if (precond == "diagonal") {
-			LOG4CPLUS_INFO(logger, "Using diagonal preconditioner.");
-		} else if (precond == "ildlt") {
-			LOG4CPLUS_INFO(logger, "Using ILDLT preconditioner.");
-		} else if (precond == "none") {
-			LOG4CPLUS_INFO(logger, "Using no preconditioner.");
+		LOG4CPLUS_INFO(logger, "Starting iterative solver.");
+		if (s->getString("lemethod") == "jacobi") {
+			if (precond != "none") {
+				LOG4CPLUS_WARN(logger, "Requested preconditioner '" << precond << "', which is unavailable for the Jacobi method. Dropping preconditioner.");
+			}
+		} else {
+			if (precond == "ilu") {
+				LOG4CPLUS_INFO(logger, "Using ILU preconditioner.");
+			} else if (precond == "diagonal") {
+				LOG4CPLUS_INFO(logger, "Using diagonal preconditioner.");
+			} else if (precond == "ildlt") {
+				LOG4CPLUS_INFO(logger, "Using ILDLT preconditioner.");
+			} else if (precond == "none") {
+				LOG4CPLUS_INFO(logger, "Using no preconditioner.");
+			}
 		}
 
+		// Now do the actual solving.
 		if (s->getString("lemethod") == "bicgstab") {
 			LOG4CPLUS_INFO(logger, "Using BiCGStab method.");
 			// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-			gmm::csr_matrix<Type>* A = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(matrix);
+			gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
 			if (precond == "ilu") {
-				gmm::bicgstab(*A, vector, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::bicgstab(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
 			} else if (precond == "diagonal") {
-				gmm::bicgstab(*A, vector, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
 			} else if (precond == "ildlt") {
-				gmm::bicgstab(*A, vector, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::bicgstab(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
 			} else if (precond == "none") {
-				gmm::bicgstab(*A, vector, b, gmm::identity_matrix(), iter);
+				gmm::bicgstab(*gmmA, x, b, gmm::identity_matrix(), iter);
 			}
 
 			// Check if the solver converged and issue a warning otherwise.
@@ -178,19 +186,19 @@ private:
 			} else {
 				LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
 			}
-			delete A;
+			delete gmmA;
 		} else if (s->getString("lemethod") == "qmr") {
 			LOG4CPLUS_INFO(logger, "Using QMR method.");
 			// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-			gmm::csr_matrix<Type>* A = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(matrix);
+			gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
 			if (precond == "ilu") {
-				gmm::qmr(*A, vector, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::qmr(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
 			} else if (precond == "diagonal") {
-				gmm::qmr(*A, vector, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::qmr(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
 			} else if (precond == "ildlt") {
-				gmm::qmr(*A, vector, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::qmr(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
 			} else if (precond == "none") {
-				gmm::qmr(*A, vector, b, gmm::identity_matrix(), iter);
+				gmm::qmr(*gmmA, x, b, gmm::identity_matrix(), iter);
 			}
 
 			// Check if the solver converged and issue a warning otherwise.
@@ -199,58 +207,62 @@ private:
 			} else {
 				LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
 			}
-			delete A;
+			delete gmmA;
 		} else if (s->getString("lemethod") == "jacobi") {
 			LOG4CPLUS_INFO(logger, "Using Jacobi method.");
-			solveLinearEquationSystemWithJacobi(matrix, vector, b);
+			uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b);
+
+			// Check if the solver converged and issue a warning otherwise.
+			if (iterations < s->get<unsigned>("maxiter")) {
+				LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
+			} else {
+				LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
+			}
 		}
 	}
 
 	/*!
-	 * Solves the linear equation system Ax=b with the given parameters
-	 * using the Jacobi Method and therefor the Jacobi Decomposition of A.
+	 * Solves the linear equation system A*x = b given by the parameters using the Jacobi method.
 	 *
-	 * @param A The matrix A specifying the coefficients of the linear equations.
-	 * @param x The vector x for which to solve the equations. The initial value of the elements of
-	 * this vector are used as the initial guess and might thus influence performance and convergence.
-	 * @param b The vector b specifying the values on the right-hand-sides of the equations.
-	 * @return The solution of the system of linear equations in form of the elements of the vector
-	 * x.
+	 * @param A The matrix specifying the coefficients of the linear equations.
+	 * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but
+	 * may be ignored.
+	 * @param b The right-hand side of the equation system.
+	 * @returns The solution vector x of the system of linear equations as the content of the parameter x.
+	 * @returns The number of iterations needed until convergence.
 	 */
-	void solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& x, std::vector<Type> const& b) const {
+	uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const {
 		// Get the settings object to customize linear solving.
 		storm::settings::Settings* s = storm::settings::instance();
 
 		double precision = s->get<double>("precision");
-		if (precision <= 0) {
-			LOG4CPLUS_ERROR(logger, "Selected precision for linear equation solving must be strictly greater than zero for Jacobi method.");
-		}
-
-		// Get a Jacobi Decomposition of the Input Matrix
-		storm::storage::JacobiDecomposition<Type>* jacobiDecomposition = matrix.getJacobiDecomposition();
+		uint_fast64_t maxIterations = s->get<unsigned>("maxiter");
+		bool relative = s->get<bool>("relative");
 
-		// Convert the Input Matrix to GMM Format so we can calculate the Residuum
-		gmm::csr_matrix<Type>* A = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(matrix);
+		// Get a Jacobi decomposition of the matrix A.
+		storm::storage::JacobiDecomposition<Type>* jacobiDecomposition = A.getJacobiDecomposition();
 
-		// Convert the Diagonal matrix to GMM format
+		// Convert the (inverted) diagonal matrix to gmm++'s format.
 		gmm::csr_matrix<Type>* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(jacobiDecomposition->getJacobiDInvReference());
-		// Convert the LU Matrix to GMM format
+		// Convert the LU matrix to gmm++'s format.
 		gmm::csr_matrix<Type>* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(jacobiDecomposition->getJacobiLUReference());
 
-
 		LOG4CPLUS_INFO(logger, "Starting iterative Jacobi Solver.");
 		
 		// x_(k + 1) = D^-1 * (b  - R * x_k)
-		// In x we keep a copy of the result for swapping in the loop (e.g. less copy-back)
+		// In x we keep a copy of the result for swapping in the loop (e.g. less copy-back).
 		std::vector<Type>* xNext = new std::vector<Type>(x.size());
 		const std::vector<Type>* xCopy = xNext;
 		std::vector<Type>* xCurrent = &x;
 
-		// Target vector for precision calculation
+		// Target vector for precision calculation.
 		std::vector<Type>* residuum = new std::vector<Type>(x.size());
 
+		// Set up additional environment variables.
 		uint_fast64_t iterationCount = 0;
-		do {
+		bool converged = false;
+
+		while (!converged && iterationCount < maxIterations) {
 			// R * x_k (xCurrent is x_k) -> xNext
 			gmm::mult(*gmmxxLU, *xCurrent, *xNext);
 			// b - R * x_k (xNext contains R * x_k) -> xNext
@@ -258,39 +270,39 @@ private:
 			// D^-1 * (b - R * x_k) -> xNext
 			gmm::mult(*gmmxxDiagonalInverted, *xNext, *xNext);
 			
-			// swap xNext with xCurrent so that the next iteration can use xCurrent again without having to copy the vector
+			// Swap xNext with xCurrent so that the next iteration can use xCurrent again without having to copy the
+			// vector.
 			std::vector<Type> *const swap = xNext;
 			xNext = xCurrent;
 			xCurrent = swap;
 
+			// Now check if the process already converged within our precision.
+			converged = storm::utility::equalModuloPrecision(*xCurrent, *xNext, precision, relative);
+
+			// Increase iteration count so we can abort if convergence is too slow.
 			++iterationCount;
-			// Precision calculation via ||A * x_k - b|| < precision
-			gmm::mult(*A, *xCurrent, *residuum);
-			gmm::add(gmm::scaled(*residuum, -1.0), b, *residuum);
-		} while (gmm::vect_norminf(*residuum) > precision);
+		}
 
-		// If the last iteration did not write to the original x
-		// we have to swith them
+		// If the last iteration did not write to the original x we have to swap the contents, because the output has to
+		// be written to the parameter x.
 		if (xCurrent == xCopy) {
-			x.swap(*xCurrent);
+			std::swap(x, *xCurrent);
 		}
 
-		// xCopy always points to the Swap-Copy of x we created
+		// As xCopy always points to the copy of x used for swapping, we can safely delete it.
 		delete xCopy;
-		// Delete the residuum vector
+
+		// Also delete the other dynamically allocated variables.
 		delete residuum;
-		// Delete the decompositions
 		delete jacobiDecomposition;
-		// and the GMM Matrices
 		delete gmmxxDiagonalInverted;
 		delete gmmxxLU;
 
-		LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterationCount << " iterations.");
+		return iterationCount;
 	}
 };
 
-} //namespace modelChecker
-
-} //namespace storm
+} // namespace modelchecker
+} // namespace storm
 
 #endif /* STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_ */
diff --git a/src/modelchecker/GmmxxMdpPrctlModelChecker.h b/src/modelchecker/GmmxxMdpPrctlModelChecker.h
index a0c2cae85..48c4c6304 100644
--- a/src/modelchecker/GmmxxMdpPrctlModelChecker.h
+++ b/src/modelchecker/GmmxxMdpPrctlModelChecker.h
@@ -8,65 +8,73 @@
 #ifndef STORM_MODELCHECKER_GMMXXMDPPRCTLMODELCHECKER_H_
 #define STORM_MODELCHECKER_GMMXXMDPPRCTLMODELCHECKER_H_
 
-#include <cmath>
-
-#include "src/models/Mdp.h"
-#include "src/modelchecker/MdpPrctlModelChecker.h"
-#include "src/utility/GraphAnalyzer.h"
-#include "src/utility/Vector.h"
-#include "src/utility/ConstTemplates.h"
-#include "src/utility/Settings.h"
+#include "src/modelchecker/SparseMdpPrctlModelChecker.h"
 #include "src/adapters/GmmxxAdapter.h"
-#include "src/exceptions/InvalidPropertyException.h"
-#include "src/storage/JacobiDecomposition.h"
 
 #include "gmm/gmm_matrix.h"
 #include "gmm/gmm_iter_solvers.h"
 
-#include "log4cplus/logger.h"
-#include "log4cplus/loggingmacros.h"
-
-extern log4cplus::Logger logger;
+#include <cmath>
 
 namespace storm {
-
-namespace modelChecker {
+namespace modelchecker {
 
 /*
- * A model checking engine that makes use of the gmm++ backend.
+ * An implementation of the SparseMdpPrctlModelChecker interface that uses gmm++ as the solving backend.
  */
 template <class Type>
-class GmmxxMdpPrctlModelChecker : public MdpPrctlModelChecker<Type> {
+class GmmxxMdpPrctlModelChecker : public SparseMdpPrctlModelChecker<Type> {
 
 public:
-	explicit GmmxxMdpPrctlModelChecker(storm::models::Mdp<Type>& mdp) : MdpPrctlModelChecker<Type>(mdp) { }
+	/*!
+	 * Constructs a GmmxxMdpPrctlModelChecker with the given model.
+	 *
+	 * @param model The MDP to be checked.
+	 */
+	explicit GmmxxMdpPrctlModelChecker(storm::models::Mdp<Type> const& model) : SparseMdpPrctlModelChecker<Type>(model) {
+		// Intentionally left empty.
+	}
 
-	virtual ~GmmxxMdpPrctlModelChecker() { }
+	/*!
+	 * Copy constructs a GmmxxDtmcPrctlModelChecker from the given model checker. In particular, this means that the newly
+	 * constructed model checker will have the model of the given model checker as its associated model.
+	 */
+	explicit GmmxxMdpPrctlModelChecker(storm::modelchecker::GmmxxMdpPrctlModelChecker<Type> const& modelchecker) : SparseMdpPrctlModelChecker<Type>(modelchecker) {
+		// Intentionally left empty.
+	}
+
+	/*!
+	 * Virtual destructor. Needs to be virtual, because this class has virtual methods.
+	 */
+	virtual ~GmmxxMdpPrctlModelChecker() {
+		// Intentionally left empty.
+	}
 
 private:
-	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type>* summand = nullptr, uint_fast64_t repetitions = 1) const {
+
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const {
 		// Get the starting row indices for the non-deterministic choices to reduce the resulting
 		// vector properly.
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
+		std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = *this->getModel().getNondeterministicChoiceIndices();
 
 
 		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(matrix);
+		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
 
 		// Create vector for result of multiplication, which is reduced to the result vector after
 		// each multiplication.
-		std::vector<Type> multiplyResult(matrix.getRowCount());
+		std::vector<Type> multiplyResult(A.getRowCount());
 
 		// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
-		for (uint_fast64_t i = 0; i < repetitions; ++i) {
-			gmm::mult(*gmmxxMatrix, vector, multiplyResult);
-			if (summand != nullptr) {
-				gmm::add(*summand, multiplyResult);
+		for (uint_fast64_t i = 0; i < n; ++i) {
+			gmm::mult(*gmmxxMatrix, x, multiplyResult);
+			if (b != nullptr) {
+				gmm::add(*b, multiplyResult);
 			}
 			if (this->minimumOperatorStack.top()) {
-				storm::utility::reduceVectorMin(multiplyResult, &vector, *nondeterministicChoiceIndices);
+				storm::utility::reduceVectorMin(multiplyResult, &x, nondeterministicChoiceIndices);
 			} else {
-				storm::utility::reduceVectorMax(multiplyResult, &vector, *nondeterministicChoiceIndices);
+				storm::utility::reduceVectorMax(multiplyResult, &x, nondeterministicChoiceIndices);
 			}
 		}
 
@@ -84,7 +92,7 @@ private:
 	 * @return The solution of the system of linear equations in form of the elements of the vector
 	 * x.
 	 */
-	void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
+	void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
 		// Get the settings object to customize solving.
 		storm::settings::Settings* s = storm::settings::instance();
 
@@ -95,18 +103,18 @@ private:
 
 
 		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(matrix);
+		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
 
 		// Set up the environment for the power method.
-		std::vector<Type> multiplyResult(matrix.getRowCount());
+		std::vector<Type> multiplyResult(A.getRowCount());
 		std::vector<Type>* currentX = &x;
 		std::vector<Type>* newX = new std::vector<Type>(x.size());
 		std::vector<Type>* swap = nullptr;
 		uint_fast64_t iterations = 0;
 		bool converged = false;
 
-		// Proceed with the iterations as long as the method did not converge or reach the
-		// user-specified maximum number of iterations.
+		// Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number
+		// of iterations.
 		while (!converged && iterations < maxIterations) {
 			// Compute x' = A*x + b.
 			gmm::mult(*gmmxxMatrix, *currentX, multiplyResult);
@@ -129,6 +137,8 @@ private:
 			++iterations;
 		}
 
+		// If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
+		// is currently stored in currentX, but x is the output vector.
 		if (iterations % 2 == 1) {
 			std::swap(x, *currentX);
 			delete currentX;
@@ -146,8 +156,7 @@ private:
 	}
 };
 
-} //namespace modelChecker
-
+} //namespace modelchecker
 } //namespace storm
 
 #endif /* STORM_MODELCHECKER_GMMXXMDPPRCTLMODELCHECKER_H_ */
diff --git a/src/modelchecker/SparseDtmcPrctlModelChecker.h b/src/modelchecker/SparseDtmcPrctlModelChecker.h
new file mode 100644
index 000000000..4bcf4708a
--- /dev/null
+++ b/src/modelchecker/SparseDtmcPrctlModelChecker.h
@@ -0,0 +1,454 @@
+/*
+ * SparseDtmcPrctlModelChecker.h
+ *
+ *  Created on: 22.10.2012
+ *      Author: Thomas Heinemann
+ */
+
+#ifndef STORM_MODELCHECKER_SPARSEDTMCPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_SPARSEDTMCPRCTLMODELCHECKER_H_
+
+#include "src/modelchecker/AbstractModelChecker.h"
+#include "src/models/Dtmc.h"
+#include "src/utility/Vector.h"
+#include "src/utility/GraphAnalyzer.h"
+
+#include <vector>
+
+namespace storm {
+namespace modelchecker {
+
+/*!
+ * @brief
+ * Interface for all model checkers that can verify PRCTL formulae over DTMCs represented as a sparse matrix.
+ */
+template<class Type>
+class SparseDtmcPrctlModelChecker : public AbstractModelChecker<Type> {
+
+public:
+	/*!
+	 * Constructs a SparseDtmcPrctlModelChecker with the given model.
+	 *
+	 * @param model The DTMC to be checked.
+	 */
+	explicit SparseDtmcPrctlModelChecker(storm::models::Dtmc<Type> const& model) : AbstractModelChecker<Type>(model) {
+		// Intentionally left empty.
+	}
+
+	/*!
+	 * Copy constructs a SparseDtmcPrctlModelChecker from the given model checker. In particular, this means that the newly
+	 * constructed model checker will have the model of the given model checker as its associated model.
+	 */
+	explicit SparseDtmcPrctlModelChecker(storm::modelchecker::SparseDtmcPrctlModelChecker<Type> const& modelChecker) : AbstractModelChecker<Type>(modelChecker) {
+		// Intentionally left empty.
+	}
+
+	/*!
+	 * Virtual destructor. Needs to be virtual, because this class has virtual methods.
+	 */
+	virtual ~SparseDtmcPrctlModelChecker() {
+		// Intentionally left empty.
+	}
+
+	/*!
+	 * Returns a constant reference to the DTMC associated with this model checker.
+	 * @returns A constant reference to the DTMC associated with this model checker.
+	 */
+	storm::models::Dtmc<Type> const& getModel() const {
+		return AbstractModelChecker<Type>::template getModel<storm::models::Dtmc<Type>>();
+	}
+
+	/*!
+	 * Checks the given formula that is a P/R operator without a bound.
+	 *
+	 * @param formula The formula to check.
+	 * @returns The set of states satisfying the formula represented by a bit vector.
+	 */
+	std::vector<Type>* checkNoBoundOperator(storm::formula::NoBoundOperator<Type> const& formula) const {
+		// Check if the operator was an optimality operator and report a warning in that case.
+		if (formula.isOptimalityOperator()) {
+			LOG4CPLUS_WARN(logger, "Formula contains min/max operator, which is not meaningful over deterministic models.");
+		}
+		return formula.getPathFormula().check(*this, false);
+	}
+
+	/*!
+	 * Checks the given formula that is a bounded-until formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
+	 */
+	virtual std::vector<Type>* checkBoundedUntil(storm::formula::BoundedUntil<Type> const& formula, bool qualitative) const {
+		// 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);
+
+		// Copy the matrix before we make any changes.
+		storm::storage::SparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionMatrix());
+
+		// Make all rows absorbing that violate both sub-formulas or satisfy the second sub-formula.
+		tmpMatrix.makeRowsAbsorbing(~(*leftStates | *rightStates) | *rightStates);
+
+		// Delete obsolete intermediates.
+		delete leftStates;
+		delete rightStates;
+
+		// Create the vector with which to multiply.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+		storm::utility::setVectorValues(result, *rightStates, storm::utility::constGetOne<Type>());
+
+		// Perform the matrix vector multiplication as often as required by the formula bound.
+		this->performMatrixVectorMultiplication(tmpMatrix, *result, nullptr, formula.getBound());
+
+		// Return result.
+		return result;
+	}
+
+	/*!
+	 * Checks the given formula that is a next formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
+	 */
+	virtual std::vector<Type>* checkNext(storm::formula::Next<Type> const& formula, bool qualitative) const {
+		// First, we need to compute the states that satisfy the child formula of the next-formula.
+		storm::storage::BitVector* nextStates = formula.getChild().check(*this);
+
+		// Create the vector with which to multiply and initialize it correctly.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+		storm::utility::setVectorValues(result, *nextStates, storm::utility::constGetOne<Type>());
+
+		// Delete obsolete intermediate.
+		delete nextStates;
+
+		// Perform one single matrix-vector multiplication.
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result);
+
+		// Return result.
+		return result;
+	}
+
+	/*!
+	 * Checks the given formula that is a bounded-eventually formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
+	 */
+	virtual std::vector<Type>* checkBoundedEventually(storm::formula::BoundedEventually<Type> const& formula, bool qualitative) const {
+		// Create equivalent temporary bounded until formula and check it.
+		storm::formula::BoundedUntil<Type> temporaryBoundedUntilFormula(new storm::formula::Ap<Type>("true"), formula.getChild().clone(), formula.getBound());
+		return this->checkBoundedUntil(temporaryBoundedUntilFormula, qualitative);
+	}
+
+	/*!
+	 * Checks the given formula that is an eventually formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
+	 */
+	virtual std::vector<Type>* checkEventually(storm::formula::Eventually<Type> const& formula, bool qualitative) const {
+		// Create equivalent temporary until formula and check it.
+		storm::formula::Until<Type> temporaryUntilFormula(new storm::formula::Ap<Type>("true"), formula.getChild().clone());
+		return this->checkUntil(temporaryUntilFormula, qualitative);
+	}
+
+	/*!
+	 * Checks the given formula that is a globally formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
+	 */
+	virtual std::vector<Type>* checkGlobally(storm::formula::Globally<Type> const& formula, bool qualitative) const {
+		// Create "equivalent" (equivalent up to negation) temporary eventually formula and check it.
+		storm::formula::Eventually<Type> temporaryEventuallyFormula(new storm::formula::Not<Type>(formula.getChild().clone()));
+		std::vector<Type>* result = this->checkEventually(temporaryEventuallyFormula, qualitative);
+
+		// Now subtract the resulting vector from the constant one vector to obtain final result.
+		storm::utility::subtractFromConstantOneVector(result);
+		return result;
+	}
+
+	/*!
+	 * Check the given formula that is an until formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
+	 */
+	virtual std::vector<Type>* checkUntil(storm::formula::Until<Type> const& formula, bool qualitative) const {
+		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
+		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
+		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
+
+		// Then, we need to identify the states which have to be taken out of the matrix, i.e.
+		// all states that have probability 0 and 1 of satisfying the until-formula.
+		storm::storage::BitVector statesWithProbability0(this->getModel().getNumberOfStates());
+		storm::storage::BitVector statesWithProbability1(this->getModel().getNumberOfStates());
+		storm::utility::GraphAnalyzer::performProb01(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
+
+		// Delete intermediate results that are obsolete now.
+		delete leftStates;
+		delete rightStates;
+
+		// Perform some logging.
+		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::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::setVectorValues<Type>(result, maybeStates, Type(0.5));
+		}
+
+		// Set values of resulting vector that are known exactly.
+		storm::utility::setVectorValues<Type>(result, statesWithProbability0, storm::utility::constGetZero<Type>());
+		storm::utility::setVectorValues<Type>(result, statesWithProbability1, storm::utility::constGetOne<Type>());
+
+		return result;
+	}
+
+	/*!
+	 * Checks the given formula that is an instantaneous reward formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bound 0. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bound 0.
+	 * @returns The reward values for the given formula for every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact values might not be computed.
+	 */
+	virtual std::vector<Type>* checkInstantaneousReward(storm::formula::InstantaneousReward<Type> const& formula, bool qualitative) const {
+		// Only compute the result if the model has a state-based reward model.
+		if (!this->getModel().hasStateRewards()) {
+			LOG4CPLUS_ERROR(logger, "Missing (state-based) reward model for formula.");
+			throw storm::exceptions::InvalidPropertyException() << "Missing (state-based) reward model for formula.";
+		}
+
+		// Initialize result to state rewards of the model.
+		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
+
+		// 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());
+
+		// Return result.
+		return result;
+	}
+
+	/*!
+	 * Check the given formula that is a cumulative reward formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bound 0. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bound 0.
+	 * @returns The reward values for the given formula for every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact values might not be computed.
+	 */
+	virtual std::vector<Type>* checkCumulativeReward(storm::formula::CumulativeReward<Type> const& formula, bool qualitative) const {
+		// Only compute the result if the model has at least one reward model.
+		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
+			LOG4CPLUS_ERROR(logger, "Missing reward model for formula.");
+			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
+		}
+
+		// Compute the reward vector to add in each step based on the available reward models.
+		std::vector<Type> totalRewardVector;
+		if (this->getModel().hasTransitionRewards()) {
+			totalRewardVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
+			if (this->getModel().hasStateRewards()) {
+				gmm::add(*this->getModel().getStateRewardVector(), totalRewardVector);
+			}
+		} else {
+			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());
+		} 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());
+
+		// Return result.
+		return result;
+	}
+
+	/*!
+	 * Checks the given formula that is a reachability reward formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bound 0. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bound 0.
+	 * @returns The reward values for the given formula for every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact values might not be computed.
+	 */
+	virtual std::vector<Type>* checkReachabilityReward(storm::formula::ReachabilityReward<Type> const& formula, bool qualitative) const {
+		// Only compute the result if the model has at least one reward model.
+		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
+			LOG4CPLUS_ERROR(logger, "Missing reward model for formula. Skipping formula");
+			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
+		}
+
+		// Determine the states for which the target predicate holds.
+		storm::storage::BitVector* targetStates = formula.getChild().check(*this);
+
+		// Determine which states have a reward of infinity by definition.
+		storm::storage::BitVector infinityStates(this->getModel().getNumberOfStates());
+		storm::storage::BitVector trueStates(this->getModel().getNumberOfStates(), true);
+		storm::utility::GraphAnalyzer::performProb1(this->getModel(), trueStates, *targetStates, &infinityStates);
+		infinityStates.complement();
+
+		// 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::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::selectVectorValues(&subStateRewards, maybeStates, *this->getModel().getStateRewardVector());
+					gmm::add(subStateRewards, b);
+				}
+			} else {
+				// If only a state-based reward model is  available, we take this vector as the
+				// right-hand side. As the state reward vector contains entries not just for the
+				// states that we still consider (i.e. maybeStates), we need to extract these values
+				// first.
+				storm::utility::selectVectorValues(&b, maybeStates, *this->getModel().getStateRewardVector());
+			}
+
+			// Now solve the resulting equation system.
+			this->solveEquationSystem(submatrix, x, b);
+
+			// Set values of resulting vector according to result.
+			storm::utility::setVectorValues<Type>(result, maybeStates, x);
+		}
+
+		// Set values of resulting vector that are known exactly.
+		storm::utility::setVectorValues(result, *targetStates, storm::utility::constGetZero<Type>());
+		storm::utility::setVectorValues(result, infinityStates, storm::utility::constGetInfinity<Type>());
+
+		// Delete temporary storages and return result.
+		delete targetStates;
+		return result;
+	}
+
+private:
+	/*!
+	 * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes x[i+1] = A*x[i] + b
+	 * until x[n], where x[0] = x.
+	 *
+	 * @param A The matrix that is to be multiplied against the vector.
+	 * @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter,
+	 * i.e. after the method returns, this vector will contain the computed values.
+	 * @param b If not null, this vector is being added to the result after each matrix-vector multiplication.
+	 * @param n Specifies the number of iterations the matrix-vector multiplication is performed.
+	 * @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector.
+	 */
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const = 0;
+
+	/*!
+	 * Solves the equation system A*x = b given by the parameters.
+	 *
+	 * @param A The matrix specifying the coefficients of the linear equations.
+	 * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but
+	 * may be ignored.
+	 * @param b The right-hand side of the equation system.
+	 * @returns The solution vector x of the system of linear equations as the content of the parameter x.
+	 */
+	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const = 0;
+};
+
+} // namespace modelchecker
+} // namespace storm
+
+#endif /* STORM_MODELCHECKER_SPARSEDTMCPRCTLMODELCHECKER_H_ */
diff --git a/src/modelchecker/MdpPrctlModelChecker.h b/src/modelchecker/SparseMdpPrctlModelChecker.h
similarity index 53%
rename from src/modelchecker/MdpPrctlModelChecker.h
rename to src/modelchecker/SparseMdpPrctlModelChecker.h
index acc75b2fb..df9d21d8c 100644
--- a/src/modelchecker/MdpPrctlModelChecker.h
+++ b/src/modelchecker/SparseMdpPrctlModelChecker.h
@@ -1,175 +1,70 @@
 /*
- * MdpPrctlModelChecker.h
+ * SparseMdpPrctlModelChecker.h
  *
  *  Created on: 15.02.2013
  *      Author: Christian Dehnert
  */
 
-#ifndef STORM_MODELCHECKER_MDPPRCTLMODELCHECKER_H_
-#define STORM_MODELCHECKER_MDPPRCTLMODELCHECKER_H_
-
-#include "src/formula/Formulas.h"
-#include "src/utility/Vector.h"
-#include "src/storage/SparseMatrix.h"
+#ifndef STORM_MODELCHECKER_SPARSEMDPPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_SPARSEMDPPRCTLMODELCHECKER_H_
 
+#include "src/modelchecker/AbstractModelChecker.h"
 #include "src/models/Mdp.h"
-#include "src/storage/BitVector.h"
-#include "src/exceptions/InvalidPropertyException.h"
 #include "src/utility/Vector.h"
-#include "src/modelchecker/AbstractModelChecker.h"
+#include "src/utility/GraphAnalyzer.h"
+
 #include <vector>
 #include <stack>
 
-#include "log4cplus/logger.h"
-#include "log4cplus/loggingmacros.h"
-
-extern log4cplus::Logger logger;
-
 namespace storm {
-
-namespace modelChecker {
+namespace modelchecker {
 
 /*!
  * @brief
- * Interface for model checker classes.
- *
- * This class provides basic functions that are the same for all subclasses, but mainly only declares
- * abstract methods that are to be implemented in concrete instances.
- *
- * @attention This class is abstract.
+ * Interface for all model checkers that can verify PRCTL formulae over MDPs represented as a sparse matrix.
  */
 template<class Type>
-class MdpPrctlModelChecker :
-	public AbstractModelChecker<Type> {
+class SparseMdpPrctlModelChecker : public AbstractModelChecker<Type> {
 
 public:
 	/*!
-	 * Constructor
+	 * Constructs a SparseMdpPrctlModelChecker with the given model.
 	 *
-	 * @param model The dtmc model which is checked.
+	 * @param model The MDP to be checked.
 	 */
-	explicit MdpPrctlModelChecker(storm::models::Mdp<Type>& model)
-		: AbstractModelChecker<Type>(model), minimumOperatorStack() {
+	explicit SparseMdpPrctlModelChecker(storm::models::Mdp<Type> const& model) : AbstractModelChecker<Type>(model), minimumOperatorStack() {
+		// Intentionally left empty.
 	}
 
 	/*!
-	 * Copy constructor
-	 *
-	 * @param modelChecker The model checker that is copied.
+	 * Copy constructs a SparseMdpPrctlModelChecker from the given model checker. In particular, this means that the newly
+	 * constructed model checker will have the model of the given model checker as its associated model.
 	 */
-	explicit MdpPrctlModelChecker(const storm::modelChecker::MdpPrctlModelChecker<Type>* modelChecker)
-		: AbstractModelChecker<Type>(modelChecker),  minimumOperatorStack() {
-
+	explicit SparseMdpPrctlModelChecker(storm::modelchecker::SparseMdpPrctlModelChecker<Type> const& modelchecker)
+		: AbstractModelChecker<Type>(modelchecker),  minimumOperatorStack() {
+		// Intentionally left empty.
 	}
 
 	/*!
-	 * Destructor
+	 * Virtual destructor. Needs to be virtual, because this class has virtual methods.
 	 */
-	virtual ~MdpPrctlModelChecker() {
+	virtual ~SparseMdpPrctlModelChecker() {
 		// Intentionally left empty.
 	}
 
 	/*!
-	 * @returns A reference to the dtmc of the model checker.
+	 * Returns a constant reference to the MDP associated with this model checker.
+	 * @returns A constant reference to the MDP associated with this model checker.
 	 */
-	storm::models::Mdp<Type>& getModel() const {
+	storm::models::Mdp<Type> const& getModel() const {
 		return AbstractModelChecker<Type>::template getModel<storm::models::Mdp<Type>>();
 	}
 
 	/*!
-	 * Sets the DTMC model which is checked
-	 * @param model
-	 */
-	void setModel(storm::models::Mdp<Type>& model) {
-		AbstractModelChecker<Type>::setModel(model);
-	}
-
-	/*!
-	 * Checks the given state formula on the DTMC and prints the result (true/false) for all initial
-	 * states.
-	 * @param stateFormula The formula to be checked.
-	 */
-	void check(const storm::formula::AbstractStateFormula<Type>& stateFormula) const {
-		std::cout << std::endl;
-		LOG4CPLUS_INFO(logger, "Model checking formula\t" << stateFormula.toString());
-		std::cout << "Model checking formula:\t" << stateFormula.toString() << std::endl;
-		storm::storage::BitVector* result = nullptr;
-		try {
-			result = stateFormula.check(*this);
-			LOG4CPLUS_INFO(logger, "Result for initial states:");
-			std::cout << "Result for initial states:" << std::endl;
-			for (auto initialState : *this->getModel().getLabeledStates("init")) {
-				LOG4CPLUS_INFO(logger, "\t" << initialState << ": " << (result->get(initialState) ? "satisfied" : "not satisfied"));
-				std::cout << "\t" << initialState << ": " << (*result)[initialState] << std::endl;
-			}
-			delete result;
-		} catch (std::exception& e) {
-			std::cout << "Error during computation: " << e.what() << "Skipping property." << std::endl;
-			if (result != nullptr) {
-				delete result;
-			}
-		}
-		std::cout << std::endl;
-		storm::utility::printSeparationLine(std::cout);
-	}
-
-	/*!
-	 * Checks the given operator (with no bound) on the DTMC and prints the result
-	 * (probability/rewards) for all initial states.
-	 * @param noBoundFormula The formula to be checked.
-	 */
-	void check(const storm::formula::NoBoundOperator<Type>& noBoundFormula) const {
-		std::cout << std::endl;
-		LOG4CPLUS_INFO(logger, "Model checking formula\t" << noBoundFormula.toString());
-		std::cout << "Model checking formula:\t" << noBoundFormula.toString() << std::endl;
-		std::vector<Type>* result = nullptr;
-		try {
-			result = noBoundFormula.check(*this);
-			LOG4CPLUS_INFO(logger, "Result for initial states:");
-			std::cout << "Result for initial states:" << std::endl;
-			for (auto initialState : *this->getModel().getLabeledStates("init")) {
-				LOG4CPLUS_INFO(logger, "\t" << initialState << ": " << (*result)[initialState]);
-				std::cout << "\t" << initialState << ": " << (*result)[initialState] << std::endl;
-			}
-			delete result;
-		} catch (std::exception& e) {
-			std::cout << "Error during computation: " << e.what() << " Skipping property." << std::endl;
-			if (result != nullptr) {
-				delete result;
-			}
-		}
-		std::cout << std::endl;
-		storm::utility::printSeparationLine(std::cout);
-	}
-
-	/*!
-	 * The check method for a formula with an AP node as root in its formula tree
+	 * Checks the given formula that is a P/R operator without a bound.
 	 *
-	 * @param formula The Ap state formula to check
-	 * @returns The set of states satisfying the formula, represented by a bit vector
-	 */
-	storm::storage::BitVector* checkAp(const storm::formula::Ap<Type>& formula) const {
-		if (formula.getAp().compare("true") == 0) {
-			return new storm::storage::BitVector(this->getModel().getNumberOfStates(), true);
-		} else if (formula.getAp().compare("false") == 0) {
-			return new storm::storage::BitVector(this->getModel().getNumberOfStates());
-		}
-
-		if (!this->getModel().hasAtomicProposition(formula.getAp())) {
-			LOG4CPLUS_ERROR(logger, "Atomic proposition '" << formula.getAp() << "' is invalid.");
-			throw storm::exceptions::InvalidPropertyException() << "Atomic proposition '" << formula.getAp() << "' is invalid.";
-			return nullptr;
-		}
-
-		return new storm::storage::BitVector(*this->getModel().getLabeledStates(formula.getAp()));
-	}
-
-	/*!
-	 * The check method for a state formula with a probabilistic operator node without bounds as root
-	 * in its formula tree
-	 *
-	 * @param formula The state formula to check
-	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 * @param formula The formula to check.
+	 * @returns The set of states satisfying the formula represented by a bit vector.
 	 */
 	std::vector<Type>* checkNoBoundOperator(const storm::formula::NoBoundOperator<Type>& formula) const {
 		// Check if the operator was an non-optimality operator and report an error in that case.
@@ -184,10 +79,15 @@ public:
 	}
 
 	/*!
-	 * The check method for a path formula with a Bounded Until operator node as root in its formula tree
+	 * Checks the given formula that is a bounded-until formula.
 	 *
-	 * @param formula The Bounded Until path formula to check
-	 * @returns for each state the probability that the path formula holds.
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
 	 */
 	virtual std::vector<Type>* checkBoundedUntil(const storm::formula::BoundedUntil<Type>& formula, bool qualitative) const {
 		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
@@ -213,10 +113,15 @@ public:
 	}
 
 	/*!
-	 * The check method for a path formula with a Next operator node as root in its formula tree
+	 * Checks the given formula that is a next formula.
 	 *
-	 * @param formula The Next path formula to check
-	 * @returns for each state the probability that the path formula holds.
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
 	 */
 	virtual std::vector<Type>* checkNext(const storm::formula::Next<Type>& formula, bool qualitative) const {
 		// First, we need to compute the states that satisfy the sub-formula of the next-formula.
@@ -236,11 +141,15 @@ public:
 	}
 
 	/*!
-	 * The check method for a path formula with a Bounded Eventually operator node as root in its
-	 * formula tree
+	 * Checks the given formula that is a bounded-eventually formula.
 	 *
-	 * @param formula The Bounded Eventually path formula to check
-	 * @returns for each state the probability that the path formula holds
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
 	 */
 	virtual std::vector<Type>* checkBoundedEventually(const storm::formula::BoundedEventually<Type>& formula, bool qualitative) const {
 		// Create equivalent temporary bounded until formula and check it.
@@ -249,10 +158,15 @@ public:
 	}
 
 	/*!
-	 * The check method for a path formula with an Eventually operator node as root in its formula tree
+	 * Checks the given formula that is an eventually formula.
 	 *
-	 * @param formula The Eventually path formula to check
-	 * @returns for each state the probability that the path formula holds
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
 	 */
 	virtual std::vector<Type>* checkEventually(const storm::formula::Eventually<Type>& formula, bool qualitative) const {
 		// Create equivalent temporary until formula and check it.
@@ -261,10 +175,15 @@ public:
 	}
 
 	/*!
-	 * The check method for a path formula with a Globally operator node as root in its formula tree
+	 * Checks the given formula that is a globally formula.
 	 *
-	 * @param formula The Globally path formula to check
-	 * @returns for each state the probability that the path formula holds
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
 	 */
 	virtual std::vector<Type>* checkGlobally(const storm::formula::Globally<Type>& formula, bool qualitative) const {
 		// Create "equivalent" temporary eventually formula and check it.
@@ -277,10 +196,15 @@ public:
 	}
 
 	/*!
-	 * The check method for a path formula with an Until operator node as root in its formula tree
+	 * Check the given formula that is an until formula.
 	 *
-	 * @param formula The Until path formula to check
-	 * @returns for each state the probability that the path formula holds.
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
 	 */
 	virtual std::vector<Type>* checkUntil(const storm::formula::Until<Type>& formula, bool qualitative) const {
 		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
@@ -314,27 +238,23 @@ 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::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
+			std::vector<uint_fast64_t> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
 
 			// Create vector for results for maybe states.
 			std::vector<Type> x(maybeStatesSetBitCount);
 
 			// Prepare the right-hand side of the equation system. For entry i this corresponds to
 			// the accumulated probability of going from state i to some 'yes' state.
-			std::vector<Type> b(submatrix->getRowCount());
-			this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, *this->getModel().getNondeterministicChoiceIndices(), statesWithProbability1, &b);
+			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);
+			this->solveEquationSystem(submatrix, x, b, subNondeterministicChoiceIndices);
 
 			// Set values of resulting vector according to result.
 			storm::utility::setVectorValues<Type>(result, maybeStates, x);
-
-			// Delete temporary matrix.
-			delete submatrix;
 		}
 
 		// Set values of resulting vector that are known exactly.
@@ -345,11 +265,15 @@ public:
 	}
 
 	/*!
-	 * The check method for a path formula with an Instantaneous Reward operator node as root in its
-	 * formula tree
+	 * Checks the given formula that is an instantaneous reward formula.
 	 *
-	 * @param formula The Instantaneous Reward formula to check
-	 * @returns for each state the reward that the instantaneous reward yields
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bound 0. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bound 0.
+	 * @returns The reward values for the given formula for every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact values might not be computed.
 	 */
 	virtual std::vector<Type>* checkInstantaneousReward(const storm::formula::InstantaneousReward<Type>& formula, bool qualitative) const {
 		// Only compute the result if the model has a state-based reward model.
@@ -368,11 +292,15 @@ public:
 	}
 
 	/*!
-	 * The check method for a path formula with a Cumulative Reward operator node as root in its
-	 * formula tree
+	 * Check the given formula that is a cumulative reward formula.
 	 *
-	 * @param formula The Cumulative Reward formula to check
-	 * @returns for each state the reward that the cumulative reward yields
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bound 0. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bound 0.
+	 * @returns The reward values for the given formula for every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact values might not be computed.
 	 */
 	virtual std::vector<Type>* checkCumulativeReward(const storm::formula::CumulativeReward<Type>& formula, bool qualitative) const {
 		// Only compute the result if the model has at least one reward model.
@@ -382,31 +310,40 @@ public:
 		}
 
 		// Compute the reward vector to add in each step based on the available reward models.
-		std::vector<Type>* totalRewardVector = nullptr;
+		std::vector<Type> totalRewardVector;
 		if (this->getModel().hasTransitionRewards()) {
 			totalRewardVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
 			if (this->getModel().hasStateRewards()) {
-				gmm::add(*this->getModel().getStateRewardVector(), *totalRewardVector);
+				gmm::add(*this->getModel().getStateRewardVector(), totalRewardVector);
 			}
 		} else {
-			totalRewardVector = new std::vector<Type>(*this->getModel().getStateRewardVector());
+			totalRewardVector = std::vector<Type>(*this->getModel().getStateRewardVector());
 		}
 
-		std::vector<Type>* result = new 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());
+		} else {
+			result = new std::vector<Type>(this->getModel().getNumberOfStates());
+		}
 
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, totalRewardVector, formula.getBound());
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, &totalRewardVector, formula.getBound());
 
 		// Delete temporary variables and return result.
-		delete totalRewardVector;
 		return result;
 	}
 
 	/*!
-	 * The check method for a path formula with a Reachability Reward operator node as root in its
-	 * formula tree
+	 * Checks the given formula that is a reachability reward formula.
 	 *
-	 * @param formula The Reachbility Reward formula to check
-	 * @returns for each state the reward that the reachability reward yields
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bound 0. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bound 0.
+	 * @returns The reward values for the given formula for every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact values might not be computed.
 	 */
 	virtual std::vector<Type>* checkReachabilityReward(const storm::formula::ReachabilityReward<Type>& formula, bool qualitative) const {
 		// Only compute the result if the model has at least one reward model.
@@ -441,25 +378,24 @@ 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::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
+			std::vector<uint_fast64_t> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
 
 			// Create vector for results for maybe states.
 			std::vector<Type> x(maybeStatesSetBitCount);
 
 			// Prepare the right-hand side of the equation system. For entry i this corresponds to
 			// the accumulated probability of going from state i to some 'yes' state.
-			std::vector<Type> b(submatrix->getRowCount());
+			std::vector<Type> b(submatrix.getRowCount());
 
 			if (this->getModel().hasTransitionRewards()) {
 				// If a transition-based reward model is available, we initialize the right-hand
 				// side to the vector resulting from summing the rows of the pointwise product
 				// of the transition probability matrix and the transition reward matrix.
-				std::vector<Type>* pointwiseProductRowSumVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
-				storm::utility::selectVectorValues(&b, maybeStates, *this->getModel().getNondeterministicChoiceIndices(), *pointwiseProductRowSumVector);
-				delete pointwiseProductRowSumVector;
+				std::vector<Type> pointwiseProductRowSumVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
+				storm::utility::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
@@ -480,11 +416,10 @@ public:
 			}
 
 			// Solve the corresponding system of equations.
-			this->solveEquationSystem(*submatrix, x, b, *subNondeterministicChoiceIndices);
+			this->solveEquationSystem(submatrix, x, b, subNondeterministicChoiceIndices);
 
 			// Set values of resulting vector according to result.
 			storm::utility::setVectorValues<Type>(result, maybeStates, x);
-			delete submatrix;
 		}
 
 		// Set values of resulting vector that are known exactly.
@@ -497,33 +432,62 @@ public:
 	}
 
 protected:
+	/*!
+	 * A stack used for storing whether we are currently computing min or max probabilities or rewards, respectively.
+	 * The topmost element is true if and only if we are currently computing minimum probabilities or rewards.
+	 */
 	mutable std::stack<bool> minimumOperatorStack;
 
 private:
-	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type>* summand = nullptr, uint_fast64_t repetitions = 1) const {
+	/*!
+	 * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes x[i+1] = A*x[i] + b
+	 * until x[n], where x[0] = x.
+	 *
+	 * @param A The matrix that is to be multiplied against the vector.
+	 * @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter,
+	 * i.e. after the method returns, this vector will contain the computed values.
+	 * @param b If not null, this vector is being added to the result after each matrix-vector multiplication.
+	 * @param n Specifies the number of iterations the matrix-vector multiplication is performed.
+	 * @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector.
+	 */
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const {
 		// Get the starting row indices for the non-deterministic choices to reduce the resulting
 		// vector properly.
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
+		std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = *this->getModel().getNondeterministicChoiceIndices();
 
 		// Create vector for result of multiplication, which is reduced to the result vector after
 		// each multiplication.
-		std::vector<Type> multiplyResult(matrix.getRowCount());
+		std::vector<Type> multiplyResult(A.getRowCount());
 
 		// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
-		for (uint_fast64_t i = 0; i < repetitions; ++i) {
-			matrix.multiplyWithVector(vector, multiplyResult);
-			if (summand != nullptr) {
-				gmm::add(*summand, multiplyResult);
+		for (uint_fast64_t i = 0; i < n; ++i) {
+			A.multiplyWithVector(x, multiplyResult);
+
+			// Add b if it is non-null.
+			if (b != nullptr) {
+				storm::utility::addVectors(multiplyResult, *b);
 			}
+
+			// Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
+			// element of the min/max operator stack.
 			if (this->minimumOperatorStack.top()) {
-				storm::utility::reduceVectorMin(multiplyResult, &vector, *nondeterministicChoiceIndices);
+				storm::utility::reduceVectorMin(multiplyResult, &x, nondeterministicChoiceIndices);
 			} else {
-				storm::utility::reduceVectorMax(multiplyResult, &vector, *nondeterministicChoiceIndices);
+				storm::utility::reduceVectorMax(multiplyResult, &x, nondeterministicChoiceIndices);
 			}
 		}
 	}
 
-	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
+	/*!
+	 * Solves the equation system A*x = b given by the parameters.
+	 *
+	 * @param A The matrix specifying the coefficients of the linear equations.
+	 * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but
+	 * may be ignored.
+	 * @param b The right-hand side of the equation system.
+	 * @returns The solution vector x of the system of linear equations as the content of the parameter x.
+	 */
+	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
 		// Get the settings object to customize solving.
 		storm::settings::Settings* s = storm::settings::instance();
 
@@ -533,7 +497,7 @@ private:
 		bool relative = s->get<bool>("relative");
 
 		// Set up the environment for the power method.
-		std::vector<Type> multiplyResult(matrix.getRowCount());
+		std::vector<Type> multiplyResult(A.getRowCount());
 		std::vector<Type>* currentX = &x;
 		std::vector<Type>* newX = new std::vector<Type>(x.size());
 		std::vector<Type>* swap = nullptr;
@@ -544,12 +508,11 @@ private:
 		// user-specified maximum number of iterations.
 		while (!converged && iterations < maxIterations) {
 			// Compute x' = A*x + b.
-			matrix.multiplyWithVector(*currentX, multiplyResult);
-			// matrix.multiplyAddAndReduceInPlace(nondeterministicChoiceIndices, *currentX, b, this->minimumOperatorStack.top());
+			A.multiplyWithVector(*currentX, multiplyResult);
+			storm::utility::addVectors(multiplyResult, b);
 
-			gmm::add(b, multiplyResult);
-
-			// Reduce the vector x' by applying min/max for all non-deterministic choices.
+			// Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
+			// element of the min/max operator stack.
 			if (this->minimumOperatorStack.top()) {
 				storm::utility::reduceVectorMin(multiplyResult, newX, nondeterministicChoiceIndices);
 			} else {
@@ -559,17 +522,15 @@ private:
 			// Determine whether the method converged.
 			converged = storm::utility::equalModuloPrecision(*currentX, *newX, precision, relative);
 
-
 			// Update environment variables.
 			swap = currentX;
 			currentX = newX;
 			newX = swap;
 			++iterations;
-
-			// *newX = *currentX,
-
 		}
 
+		// If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
+		// is currently stored in currentX, but x is the output vector.
 		if (iterations % 2 == 1) {
 			std::swap(x, *currentX);
 			delete currentX;
@@ -585,25 +546,41 @@ private:
 		}
 	}
 
-	std::shared_ptr<std::vector<uint_fast64_t>> computeNondeterministicChoiceIndicesForConstraint(storm::storage::BitVector constraint) const {
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
-		std::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices(new std::vector<uint_fast64_t>(constraint.getNumberOfSetBits() + 1));
+	/*!
+	 * Computes the nondeterministic choice indices vector resulting from reducing the full system to the states given
+	 * by the parameter constraint.
+	 *
+	 * @param constraint A bit vector specifying which states are kept.
+	 * @returns A vector of the nondeterministic choice indices of the subsystem induced by the given constraint.
+	 */
+	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();
+
+		// Reserve the known amount of slots for the resulting vector.
+		std::vector<uint_fast64_t> subNondeterministicChoiceIndices(constraint.getNumberOfSetBits() + 1);
 		uint_fast64_t currentRowCount = 0;
 		uint_fast64_t currentIndexCount = 1;
-		(*subNondeterministicChoiceIndices)[0] = 0;
+
+		// Set the first element as this will clearly begin at offset 0.
+		subNondeterministicChoiceIndices[0] = 0;
+
+		// Loop over all states that need to be kept and copy the relative indices of the nondeterministic choices over
+		// to the resulting vector.
 		for (auto index : constraint) {
-			(*subNondeterministicChoiceIndices)[currentIndexCount] = currentRowCount + (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index];
-			currentRowCount += (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index];
+			subNondeterministicChoiceIndices[currentIndexCount] = currentRowCount + nondeterministicChoiceIndices[index + 1] - nondeterministicChoiceIndices[index];
+			currentRowCount += nondeterministicChoiceIndices[index + 1] - nondeterministicChoiceIndices[index];
 			++currentIndexCount;
 		}
-		(*subNondeterministicChoiceIndices)[constraint.getNumberOfSetBits()] = currentRowCount;
+
+		// Put a sentinel element at the end.
+		subNondeterministicChoiceIndices[constraint.getNumberOfSetBits()] = currentRowCount;
 
 		return subNondeterministicChoiceIndices;
 	}
 };
 
-} //namespace modelChecker
-
-} //namespace storm
+} // namespace modelchecker
+} // namespace storm
 
-#endif /* STORM_MODELCHECKER_DTMCPRCTLMODELCHECKER_H_ */
+#endif /* STORM_MODELCHECKER_SPARSEMDPPRCTLMODELCHECKER_H_ */
diff --git a/src/modelchecker/TopologicalValueIterationMdpPrctlModelChecker.h b/src/modelchecker/TopologicalValueIterationMdpPrctlModelChecker.h
index 23b1ac6fc..248749624 100644
--- a/src/modelchecker/TopologicalValueIterationMdpPrctlModelChecker.h
+++ b/src/modelchecker/TopologicalValueIterationMdpPrctlModelChecker.h
@@ -8,39 +8,42 @@
 #ifndef STORM_MODELCHECKER_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_H_
 #define STORM_MODELCHECKER_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_H_
 
-#include <cmath>
-
-#include "src/models/Mdp.h"
-#include "src/modelchecker/MdpPrctlModelChecker.h"
-#include "src/utility/GraphAnalyzer.h"
-#include "src/utility/Vector.h"
-#include "src/utility/ConstTemplates.h"
-#include "src/utility/Settings.h"
-#include "src/adapters/GmmxxAdapter.h"
+#include "src/modelchecker/SparseMdpPrctlModelChecker.h"
 #include "src/exceptions/InvalidPropertyException.h"
-#include "src/storage/JacobiDecomposition.h"
-
-#include "gmm/gmm_matrix.h"
-#include "gmm/gmm_iter_solvers.h"
-
-#include "log4cplus/logger.h"
-#include "log4cplus/loggingmacros.h"
-
-extern log4cplus::Logger logger;
+#include <cmath>
 
 namespace storm {
-
-namespace modelChecker {
+namespace modelchecker {
 
 /*
- * A model checking engine that makes use of the gmm++ backend.
+ * An implementation of the SparseMdpPrctlModelChecker interface that uses topoligical value iteration for solving
+ * equation systems.
  */
 template <class Type>
-class TopologicalValueIterationMdpPrctlModelChecker : public MdpPrctlModelChecker<Type> {
+class TopologicalValueIterationMdpPrctlModelChecker : public SparseMdpPrctlModelChecker<Type> {
 
 public:
-	explicit TopologicalValueIterationMdpPrctlModelChecker(storm::models::Mdp<Type>& mdp) : MdpPrctlModelChecker<Type>(mdp) { }
+	/*!
+	 * Constructs a SparseMdpPrctlModelChecker with the given model.
+	 *
+	 * @param model The MDP to be checked.
+	 */
+	explicit TopologicalValueIterationMdpPrctlModelChecker(storm::models::Mdp<Type> const& model) : SparseMdpPrctlModelChecker<Type>(model) {
+		// Intentionally left empty.
+	}
+
+	/*!
+	 * Copy constructs a SparseMdpPrctlModelChecker from the given model checker. In particular, this means that the newly
+	 * constructed model checker will have the model of the given model checker as its associated model.
+	 */
+	explicit TopologicalValueIterationMdpPrctlModelChecker(storm::modelchecker::TopologicalValueIterationMdpPrctlModelChecker<Type> const& modelchecker)
+		: SparseMdpPrctlModelChecker<Type>(modelchecker),  minimumOperatorStack() {
+		// Intentionally left empty.
+	}
 
+	/*!
+	 * Virtual destructor. Needs to be virtual, because this class has virtual methods.
+	 */
 	virtual ~TopologicalValueIterationMdpPrctlModelChecker() { }
 
 private:
@@ -63,10 +66,10 @@ private:
 		unsigned maxIterations = s->get<unsigned>("maxiter");
 		bool relative = s->get<bool>("relative");
 
+		// Now, we need to determine the SCCs of the MDP and a topological sort.
 		std::vector<std::vector<uint_fast64_t>> stronglyConnectedComponents;
 		storm::models::GraphTransitions<Type> stronglyConnectedComponentsDependencyGraph;
 		storm::utility::GraphAnalyzer::performSccDecomposition(matrix, nondeterministicChoiceIndices, stronglyConnectedComponents, stronglyConnectedComponentsDependencyGraph);
-
 		std::vector<uint_fast64_t> topologicalSort;
 		storm::utility::GraphAnalyzer::getTopologicalSort(stronglyConnectedComponentsDependencyGraph, topologicalSort);
 
@@ -80,9 +83,12 @@ private:
 		uint_fast64_t globalIterations = 0;
 		bool converged = true;
 
+		// Iterate over all SCCs of the MDP as specified by the topological sort. This guarantees that an SCC is only
+		// solved after all SCCs it depends on have been solved.
 		for (auto sccIndexIt = topologicalSort.begin(); sccIndexIt != topologicalSort.end() && converged; ++sccIndexIt) {
 			std::vector<uint_fast64_t> const& scc = stronglyConnectedComponents[*sccIndexIt];
 
+			// For the current SCC, we need to perform value iteration until convergence.
 			localIterations = 0;
 			converged = false;
 			while (!converged && localIterations < maxIterations) {
@@ -98,6 +104,8 @@ private:
 				}
 
 				// Determine whether the method converged.
+				// TODO: It seems that the equalModuloPrecision call that compares all values should have a higher
+				// running time. In fact, it is faster. This has to be investigated.
 				// converged = storm::utility::equalModuloPrecision(*currentX, *newX, scc, precision, relative);
 				converged = storm::utility::equalModuloPrecision(*currentX, *newX, precision, relative);
 
@@ -109,13 +117,16 @@ private:
 				++globalIterations;
 			}
 
-			std::cout << "converged locally for scc of size " << scc.size() << std::endl;
-
+			// As the "number of iterations" of the full method is the maximum of the local iterations, we need to keep
+			// track of the maximum.
 			if (localIterations > currentMaxLocalIterations) {
 				currentMaxLocalIterations = localIterations;
 			}
 		}
 
+		// If we performed an odd number of global iterations, we need to swap the x and currentX, because the newest
+		// result is currently stored in currentX, but x is the output vector.
+		// TODO: Check whether this is correct or should be put into the for-loop over SCCs.
 		if (globalIterations % 2 == 1) {
 			std::swap(x, *currentX);
 			delete currentX;
@@ -132,8 +143,7 @@ private:
 	}
 };
 
-} //namespace modelChecker
-
-} //namespace storm
+} // namespace modelchecker
+} // namespace storm
 
 #endif /* STORM_MODELCHECKER_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_H_ */
diff --git a/src/models/Dtmc.h b/src/models/Dtmc.h
index eefcc2ab5..41e6bae8a 100644
--- a/src/models/Dtmc.h
+++ b/src/models/Dtmc.h
@@ -50,7 +50,7 @@ public:
 			throw storm::exceptions::InvalidArgumentException() << "Probability matrix is invalid.";
 		}
 		if (this->getTransitionRewardMatrix() != nullptr) {
-			if (!this->getTransitionRewardMatrix()->isSubmatrixOf(*(this->getTransitionMatrix()))) {
+			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.";
 			}
@@ -62,7 +62,7 @@ public:
 	 * Copy Constructor. Performs a deep copy of the given DTMC.
 	 * @param dtmc A reference to the DTMC that is to be copied.
 	 */
-	Dtmc(const Dtmc<T> &dtmc) : AbstractDeterministicModel<T>(dtmc) {
+	Dtmc(Dtmc<T> const& dtmc) : AbstractDeterministicModel<T>(dtmc) {
 		// Intentionally left empty.
 	}
 
@@ -95,7 +95,7 @@ private:
 			return false;
 		}
 
-		for (uint_fast64_t row = 0; row < this->getTransitionMatrix()->getRowCount(); 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");
diff --git a/src/models/GraphTransitions.h b/src/models/GraphTransitions.h
index ec8bb2bfd..feed1b952 100644
--- a/src/models/GraphTransitions.h
+++ b/src/models/GraphTransitions.h
@@ -66,29 +66,6 @@ public:
 		// Intentionally left empty.
 	}
 
-	GraphTransitions(GraphTransitions<T>&& other) {
-		this->numberOfStates = other.numberOfStates;
-		this->numberOfTransitions = other.numberOfTransitions;
-		std::swap(successorList, other.successorList);
-		std::swap(stateIndications, other.stateIndications);
-	}
-
-	GraphTransitions<T>& operator=(GraphTransitions<T>&& other) {
-		this->numberOfStates = other.numberOfStates;
-		this->numberOfTransitions = other.numberOfTransitions;
-		std::swap(successorList, other.successorList);
-		std::swap(stateIndications, other.stateIndications);
-		return *this;
-	}
-
-	//! Destructor
-	/*!
-	 * Destructor. Frees the internal storage.
-	 */
-	~GraphTransitions() {
-		// Intentionally left empty.
-	}
-
 	/*!
 	 * Retrieves the size of the internal representation of the graph transitions in memory.
 	 * @return the size of the internal representation of the graph transitions in memory
@@ -181,19 +158,24 @@ private:
 	void initializeForward(storm::storage::SparseMatrix<T> const& transitionMatrix) {
 		// First, we copy the index list from the sparse matrix as this will
 		// stay the same.
-		std::copy(transitionMatrix.getRowIndications().begin(), transitionMatrix.getRowIndications().end(), this->stateIndications.begin());
+		std::copy(transitionMatrix.constColumnIteratorBegin(), transitionMatrix.constColumnIteratorEnd(), this->stateIndications.begin());
 
 		// Now we can iterate over all rows of the transition matrix and record the target state.
 		for (uint_fast64_t i = 0, currentNonZeroElement = 0; i < numberOfStates; i++) {
-			for (auto rowIt = transitionMatrix.beginConstColumnIterator(i); rowIt != transitionMatrix.endConstColumnIterator(i); ++rowIt) {
+			for (auto rowIt = transitionMatrix.constColumnIteratorBegin(i); rowIt != transitionMatrix.constColumnIteratorEnd(i); ++rowIt) {
 				this->successorList[currentNonZeroElement++] = *rowIt;
 			}
 		}
 	}
 
 	void initializeForward(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) {
+		// We can directly copy the starting indices from the transition matrix as we do not
+		// eliminate duplicate transitions and therefore will have as many non-zero entries as this
+		// matrix.
+		typename storm::storage::SparseMatrix<T>::ConstRowsIterator rowsIt(transitionMatrix);
 		for (uint_fast64_t i = 0; i < numberOfStates; ++i) {
-			this->stateIndications[i] = transitionMatrix.getRowIndications().at(nondeterministicChoiceIndices[i]);
+			rowsIt.moveToRow(nondeterministicChoiceIndices[i]);
+			this->stateIndications[i] = rowsIt.index();
 		}
 		this->stateIndications[numberOfStates] = numberOfTransitions;
 
@@ -201,7 +183,7 @@ private:
 		// the target state.
 		for (uint_fast64_t i = 0, currentNonZeroElement = 0; i < numberOfStates; i++) {
 			for (uint_fast64_t j = nondeterministicChoiceIndices[i]; j < nondeterministicChoiceIndices[i + 1]; ++j) {
-				for (auto rowIt = transitionMatrix.beginConstColumnIterator(j); rowIt != transitionMatrix.endConstColumnIterator(j); ++rowIt) {
+				for (auto rowIt = transitionMatrix.constColumnIteratorBegin(j); rowIt != transitionMatrix.constColumnIteratorEnd(j); ++rowIt) {
 					this->successorList[currentNonZeroElement++] = *rowIt;
 				}
 			}
@@ -216,7 +198,7 @@ private:
 	void initializeBackward(storm::storage::SparseMatrix<T> const& transitionMatrix) {
 		// First, we need to count how many backward transitions each state has.
 		for (uint_fast64_t i = 0; i < numberOfStates; ++i) {
-			for (auto rowIt = transitionMatrix.beginConstColumnIterator(i); rowIt != transitionMatrix.endConstColumnIterator(i); ++rowIt) {
+			for (auto rowIt = transitionMatrix.constColumnIteratorBegin(i); rowIt != transitionMatrix.constColumnIteratorEnd(i); ++rowIt) {
 				this->stateIndications[*rowIt + 1]++;
 			}
 		}
@@ -238,7 +220,7 @@ private:
 		// Now we are ready to actually fill in the list of predecessors for
 		// every state. Again, we start by considering all but the last row.
 		for (uint_fast64_t i = 0; i < numberOfStates; ++i) {
-			for (auto rowIt = transitionMatrix.beginConstColumnIterator(i); rowIt != transitionMatrix.endConstColumnIterator(i); ++rowIt) {
+			for (auto rowIt = transitionMatrix.constColumnIteratorBegin(i); rowIt != transitionMatrix.constColumnIteratorEnd(i); ++rowIt) {
 				this->successorList[nextIndices[*rowIt]++] = i;
 			}
 		}
@@ -248,7 +230,7 @@ private:
 		// First, we need to count how many backward transitions each state has.
 		for (uint_fast64_t i = 0; i < numberOfStates; ++i) {
 			for (uint_fast64_t j = nondeterministicChoiceIndices[i]; j < nondeterministicChoiceIndices[i + 1]; ++j) {
-				for (auto rowIt = transitionMatrix.beginConstColumnIterator(j); rowIt != transitionMatrix.endConstColumnIterator(j); ++rowIt) {
+				for (auto rowIt = transitionMatrix.constColumnIteratorBegin(j); rowIt != transitionMatrix.constColumnIteratorEnd(j); ++rowIt) {
 					this->stateIndications[*rowIt + 1]++;
 				}
 			}
@@ -272,7 +254,7 @@ private:
 		// every state. Again, we start by considering all but the last row.
 		for (uint_fast64_t i = 0; i < numberOfStates; i++) {
 			for (uint_fast64_t j = nondeterministicChoiceIndices[i]; j < nondeterministicChoiceIndices[i + 1]; ++j) {
-				for (auto rowIt = transitionMatrix.beginConstColumnIterator(j); rowIt != transitionMatrix.endConstColumnIterator(j); ++rowIt) {
+				for (auto rowIt = transitionMatrix.constColumnIteratorBegin(j); rowIt != transitionMatrix.constColumnIteratorEnd(j); ++rowIt) {
 					this->successorList[nextIndices[*rowIt]++] = i;
 				}
 			}
diff --git a/src/parser/PrctlFileParser.cpp b/src/parser/PrctlFileParser.cpp
index 9efefe1a4..90721ee8f 100644
--- a/src/parser/PrctlFileParser.cpp
+++ b/src/parser/PrctlFileParser.cpp
@@ -18,16 +18,16 @@ namespace parser {
 
 PrctlFileParser::PrctlFileParser(std::string filename, storm::models::Dtmc<double>& dtmc, enum libraries library) {
 
-	storm::modelChecker::DtmcPrctlModelChecker<double>* modelChecker = nullptr;
+	storm::modelchecker::SparseDtmcPrctlModelChecker<double>* modelChecker = nullptr;
 	switch(library) {
 	//case EIGEN:
 		//Eigen Model Checker is not completely implemented at the moment, thus using Eigen is not possible...
 		//Current behaviour: Fall back to GMMXX...
-		//modelChecker = new storm::modelChecker::EigenDtmcPrctlModelChecker<double>(dtmc);
+		//modelChecker = new storm::modelchecker::EigenDtmcPrctlModelChecker<double>(dtmc);
 	//	break;
 	case GMMXX:
 	default:					//Note: GMMXX is default, hence default branches here, too.
-		modelChecker = new storm::modelChecker::GmmxxDtmcPrctlModelChecker<double>(dtmc);
+		modelChecker = new storm::modelchecker::GmmxxDtmcPrctlModelChecker<double>(dtmc);
 		break;
 	}
 	check(filename, modelChecker);
@@ -35,7 +35,7 @@ PrctlFileParser::PrctlFileParser(std::string filename, storm::models::Dtmc<doubl
 }
 
 PrctlFileParser::PrctlFileParser(std::string filename, storm::models::Mdp<double>& mdp, enum libraries library) {
-	storm::modelChecker::MdpPrctlModelChecker<double>* modelChecker = nullptr;
+	storm::modelchecker::SparseMdpPrctlModelChecker<double>* modelChecker = nullptr;
 	switch(library) {
 	//case EIGEN:
 		//Eigen MDP Model Checker is not implemented yet
@@ -43,7 +43,7 @@ PrctlFileParser::PrctlFileParser(std::string filename, storm::models::Mdp<double
 	//	break;
 	case GMMXX:
 	default:					//Note: GMMXX is default, hence default branches here, too.
-		modelChecker = new storm::modelChecker::GmmxxMdpPrctlModelChecker<double>(mdp);
+		modelChecker = new storm::modelchecker::GmmxxMdpPrctlModelChecker<double>(mdp);
 		break;
 	}
 	check(filename, modelChecker);
@@ -54,7 +54,7 @@ PrctlFileParser::~PrctlFileParser() {
 	//intentionally left empty
 }
 
-void PrctlFileParser::check(std::string filename, storm::modelChecker::AbstractModelChecker<double>* modelChecker) {
+void PrctlFileParser::check(std::string filename, storm::modelchecker::AbstractModelChecker<double>* modelChecker) {
 	// Open file
 	std::ifstream inputFileStream(filename, std::ios::in);
 
diff --git a/src/parser/PrctlFileParser.h b/src/parser/PrctlFileParser.h
index e0d64c499..5e8645c86 100644
--- a/src/parser/PrctlFileParser.h
+++ b/src/parser/PrctlFileParser.h
@@ -57,7 +57,7 @@ protected:
 	 * @param filename     The name of the file to parse
 	 * @param modelChecker The model checker that checks the formula (has to know its model!)
 	 */
-	void check(std::string filename, storm::modelChecker::AbstractModelChecker<double>* modelChecker);
+	void check(std::string filename, storm::modelchecker::AbstractModelChecker<double>* modelChecker);
 
 	/*!
 	 * Destructor.
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index 5259a8052..6dfeabf2a 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -38,8 +38,8 @@ namespace storage {
 
 /*!
  * A sparse matrix class with a constant number of non-zero entries.
- * NOTE: Addressing *is* zero-based, so the valid range for getValue and addNextValue is 0..(rows - 1)
- * where rows is the first argument to the constructor.
+ * NOTE: Addressing *is* zero-based, so the valid range for getValue and addNextValue is 0..(rows - 1) where rows is the
+ * first argument to the constructor.
  */
 template<class T>
 class SparseMatrix {
@@ -52,58 +52,148 @@ public:
 	friend class storm::adapters::StormAdapter;
 
 	/*!
-	 * If we only want to iterate over the columns of the non-zero entries of
+	 * If we only want to iterate over the columns or values of the non-zero entries of
 	 * a row, we can simply iterate over the array (part) itself.
 	 */
-	typedef const uint_fast64_t * const constIndexIterator;
+	typedef const uint_fast64_t* const ConstIndexIterator;
+	typedef const T* const ConstValueIterator;
 
 	/*!
-	 *	Iterator type if we want to iterate over elements.
+	 * Iterator class that is able to iterate over the non-zero elements of a matrix and return the position
+	 * (row, column) in addition to the value itself. This is a const iterator in the sense that it is not possible to
+	 * modify the matrix with it.
 	 */
-	typedef const T* const constIterator;
-
-	class constRowsIterator {
+	class ConstRowsIterator {
 	public:
-		constRowsIterator(SparseMatrix<T> const& matrix) : matrix(matrix), offset(0) {
+		/*!
+		 * Constructs an iterator over the elements of the given matrix.
+		 *
+		 * @param matrix The matrix on which this iterator operates.
+		 */
+		ConstRowsIterator(SparseMatrix<T> const& matrix, uint_fast64_t rowIndex = 0) : matrix(matrix), posIndex(0), rowIndex(0) {
 			// Intentionally left empty.
 		}
 
-		constRowsIterator& operator++() {
-			++offset;
+		/*!
+		 * Moves the iterator the next non-zero element of the matrix.
+		 *
+		 * @returns A reference to itself.
+		 */
+		ConstRowsIterator& operator++() {
+			++posIndex;
+			if (posIndex >= matrix.rowIndications[rowIndex + 1]) {
+				++rowIndex;
+			}
+			return *this;
+		}
+		
+		/*!
+		 * Dereferencing operator for this iterator. Actually returns a reference to itself. This is
+		 * needed, because the range-based for-loop in C++11 dereferences the iterator automatically.
+		 *
+		 * @returns A reference to itself.
+		 */
+		ConstRowsIterator& operator*() {
 			return *this;
 		}
 
-		bool operator==(constRowsIterator const& other) const {
-			return offset == other.offset;
+		/*!
+		 * Comparison operator that compares the current iterator with the given one in terms of
+         * the indices they point to. Note that this does not check whether the iterators are
+         * interpreted over the same matrix.
+		 *
+		 * @return True iff the given iterator points to the same index as the current iterator.
+		 */
+		bool operator==(ConstRowsIterator const& other) const {
+			return this->posIndex == other.posIndex;
 		}
 
-		bool operator!=(uint_fast64_t offset) const {
-			return this->offset != offset;
+		/*!
+		 * Comparison operator that compares the current iterator with the given one in terms of
+         * the indices they point to. Note that this does not check whether the iterators are
+         * interpreted over the same matrix.
+		 *
+		 * @return True iff the given iterator points to a differnent index as the current iterator.
+		 */
+		bool operator!=(ConstRowsIterator const& other) const {
+			return this->posIndex != other.posIndex;
 		}
 
-		uint_fast64_t column() const {
-			return matrix.columnIndications[offset];
+        /*!
+         * Retrieves the row that is associated with the current non-zero element this iterator
+         * points to.
+		 *
+		 * @returns The row of the current non-zero element this iterator points to.
+         */
+		uint_fast64_t row() const {
+			return this->rowIndex;
 		}
 
+        /*!
+         * Retrieves the column that is associated with the current non-zero element this iterator
+         * points to.
+		 *
+		 * @returns The column of the current non-zero element this iterator points to.
+         */
+		uint_fast64_t column() const {
+			return matrix.columnIndications[posIndex];
+		}
+		
+		/*!
+		 * Retrieves the internal index of this iterator. This index corresponds to the position of
+		 * the current element in the vector of non-zero values of the matrix.
+		 */
+		uint_fast64_t index() const {
+			return this->posIndex;
+		}
+		
+        /*!
+         * Retrieves the value of the current non-zero element this iterator points to.
+		 *
+		 * @returns The value of the current non-zero element this iterator points to.
+         */
 		T const& value() const {
-			return matrix.valueStorage[offset];
+			return matrix.valueStorage[posIndex];
+		}
+
+        /*!
+         * Moves the iterator to the beginning of the given row.
+         *
+         * @param row The row this iterator is to be moved to.
+         */
+		void moveToRow(uint_fast64_t row) {
+			this->rowIndex = row;
+			this->posIndex = matrix.rowIndications[row];
 		}
 
-		void setOffset(uint_fast64_t offset) {
-			this->offset = offset;
+        /*!
+         * Moves the iterator to the beginning of the next row.
+         */
+		void moveToNextRow() {
+			moveToRow(rowIndex + 1);
 		}
 
 	private:
+        // A constant reference to the matrix this iterator is associated with.
 		SparseMatrix<T> const& matrix;
-		uint_fast64_t offset;
+        
+        // The current index in the list of all non-zero elements of the matrix this iterator points to.
+		uint_fast64_t posIndex;
+        
+        // The row of the element this iterator currently points to.
+		uint_fast64_t rowIndex;
 	};
 
-	friend class constRowsIterator;
+    // Declare the iterator as a friend class to grant access to private data members of the matrix.
+	friend class ConstRowsIterator;
 
 	/*!
-	 * An enum representing the internal state of the Matrix.
-	 * After creating the Matrix using the Constructor, the Object is in state UnInitialized. After calling initialize(), that state changes to Initialized and after all entries have been entered and finalize() has been called, to ReadReady.
+	 * An enum representing the internal state of the Matrix
+	 * After creating the Matrix using the Constructor, the Object is in state UnInitialized. After
+	 * calling initialize(), that state changes to Initialized and after all entries have been
+	 * entered and finalize() has been called, to ReadReady.
 	 * Should a critical error occur in any of the former functions, the state will change to Error.
+	 *
 	 * @see getState()
 	 * @see isReadReady()
 	 * @see isInitialized()
@@ -113,81 +203,34 @@ public:
 		Error = -1, UnInitialized = 0, Initialized = 1, ReadReady = 2
 	};
 
-	//! Constructor
 	/*!
 	 * Constructs a sparse matrix object with the given number of rows.
+     *
 	 * @param rows The number of rows of the matrix
 	 */
-	SparseMatrix(uint_fast64_t rows, uint_fast64_t cols)
-			: rowCount(rows), colCount(cols), nonZeroEntryCount(0),
-			  internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { }
-
-	/* Sadly, Delegate Constructors are not yet available with MSVC2012 */
-	//! Constructor
-	/*!
-	 * Constructs a square sparse matrix object with the given number rows
-	 * @param size The number of rows and cols in the matrix
-	 */ /*
-	SparseMatrix(uint_fast64_t size) : SparseMatrix(size, size) { }
-	*/
+	SparseMatrix(uint_fast64_t rows, uint_fast64_t cols) : rowCount(rows), colCount(cols),
+		nonZeroEntryCount(0), internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) {
+    
+    }
 
-	//! Constructor
 	/*!
-	 * Constructs a square sparse matrix object with the given number rows
-	 * @param size The number of rows and cols in the matrix
+	 * Constructs a square sparse matrix object with the given number rows.
+	 *
+	 * @param size The number of rows and columns of the matrix.
 	 */
 	SparseMatrix(uint_fast64_t size)
 			: rowCount(size), colCount(size), nonZeroEntryCount(0),
-			  internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { }
-
-	//! Copy Constructor
-	/*!
-	 * Copy Constructor. Performs a deep copy of the given sparse matrix.
-	 * @param ssm A reference to the matrix to be copied.
-	 */
-	SparseMatrix(const SparseMatrix<T> &ssm)
-			: rowCount(ssm.rowCount), colCount(ssm.colCount), nonZeroEntryCount(ssm.nonZeroEntryCount),
-			  internalStatus(ssm.internalStatus), currentSize(ssm.currentSize), lastRow(ssm.lastRow) {
-		LOG4CPLUS_WARN(logger, "Invoking copy constructor.");
-		// Check whether copying the matrix is safe.
-		if (ssm.hasError()) {
-			LOG4CPLUS_ERROR(logger, "Trying to copy sparse matrix in error state.");
-			throw storm::exceptions::InvalidArgumentException("Trying to copy sparse matrix in error state.");
-		} else {
-			// Try to prepare the internal storage and throw an error in case
-			// of a failure.
-			if (!prepareInternalStorage()) {
-				LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage.");
-				throw std::bad_alloc();
-			} else {
-				std::copy(ssm.valueStorage.begin(), ssm.valueStorage.end(), valueStorage.begin());
-
-				// The elements that are not of the value type but rather the
-				// index type may be copied directly.
-				std::copy(ssm.columnIndications.begin(), ssm.columnIndications.end(), columnIndications.begin());
-				std::copy(ssm.rowIndications.begin(), ssm.rowIndications.end(), rowIndications.begin());
-			}
-		}
-	}
-
-	//! Destructor
-	/*!
-	 * Destructor. Performs deletion of the reserved storage arrays.
-	 */
-	~SparseMatrix() {
-		setState(MatrixStatus::UnInitialized);
-		valueStorage.resize(0);
-		columnIndications.resize(0);
-		rowIndications.resize(0);
+			  internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) {
+			  
 	}
 
 	/*!
 	 * Initializes the sparse matrix with the given number of non-zero entries
 	 * and prepares it for use with addNextValue() and finalize().
-	 * NOTE: Calling this method before any other member function is mandatory.
+	 * Note: Calling this method before any other member function is mandatory.
 	 * This version is to be used together with addNextValue().
-	 * @param nonZeroEntries The number of non-zero entries that are not on the
-	 * diagonal.
+	 *
+	 * @param nonZeroEntries The number of non-zero entries this matrix is going to hold.
 	 */
 	void initialize(uint_fast64_t nonZeroEntries) {
 		// Check whether initializing the matrix is safe.
@@ -221,17 +264,18 @@ public:
 	/*!
 	 * Sets the matrix element at the given row and column to the given value.
 	 * NOTE: This is a linear setter. It must be called consecutively for each element,
-	 * row by row *and* column by column. Only diagonal entries may be set at any time.
+	 * row by row *and* column by column.
+	 *
 	 * @param row The row in which the matrix element is to be set.
 	 * @param col The column in which the matrix element is to be set.
 	 * @param value The value that is to be set.
 	 */
-	void addNextValue(const uint_fast64_t row, const uint_fast64_t col,	const T& value) {
+	void addNextValue(const uint_fast64_t row, const uint_fast64_t col,	T const& value) {
 		// Check whether the given row and column positions are valid and throw
 		// error otherwise.
 		if ((row > rowCount) || (col > colCount)) {
 			triggerErrorState();
-			LOG4CPLUS_ERROR(logger, "Trying to add a value at illegal position (" << row << ", " << col << ").");
+			LOG4CPLUS_ERROR(logger, "Trying to add a value at illegal position (" << row << ", " << col << ") in matrix of size (" << rowCount << ", " << colCount << ").");
 			throw storm::exceptions::OutOfRangeException("Trying to add a value at illegal position.");
 		}
 
@@ -287,12 +331,13 @@ public:
 
 	/*!
 	 * Gets the matrix element at the given row and column to the given value.
-	 * NOTE: This function does not check the internal status for errors for performance reasons.
+	 * Note: This function does not check the internal status for errors for performance reasons.
+	 *
 	 * @param row The row in which the element is to be read.
 	 * @param col The column in which the element is to be read.
 	 * @param target A pointer to the memory location where the read content is
 	 * to be put.
-	 * @return True iff the value is set in the matrix, false otherwise.
+	 * @returns True iff the value is set in the matrix, false otherwise.
 	 * On false, 0 will be written to *target.
 	 */
 	inline bool getValue(uint_fast64_t row, uint_fast64_t col, T* const target) const {
@@ -314,6 +359,7 @@ public:
 				*target = valueStorage[rowStart];
 				return true;
 			}
+			
 			// If the column of the current element is already larger than the
 			// requested column, the requested element cannot be contained
 			// in the matrix and we may therefore stop searching.
@@ -329,15 +375,16 @@ public:
 	}
 
 	/*!
-	 * Gets the matrix element at the given row and column to the given value.
-	 * NOTE: This function does not check the internal status for errors for performance reasons.
-	 * WARNING: It is possible to modify through this function. Usage only valid
-	 * for elements EXISTING in the sparse matrix! If the requested value does not exist,
+	 * Gets the matrix element at the given row and column in the form of a reference to it.
+	 * Note: This function does not check the internal status for errors for performance reasons.
+	 * Warning: It is possible to modify the matrix through this function. It may only be used
+	 * for elements that exist in the sparse matrix. If the value at the requested position does not exist,
 	 * an exception will be thrown.
+	 
 	 * @param row The row in which the element is to be read.
 	 * @param col The column in which the element is to be read.
 	 *
-	 * @return A reference to the value
+	 * @return A reference to the value at the given position.
 	 */
 	inline T& getValue(uint_fast64_t row, uint_fast64_t col) {
 		// Check for illegal access indices.
@@ -351,11 +398,11 @@ public:
 		uint_fast64_t rowStart = rowIndications[row];
 		uint_fast64_t rowEnd = rowIndications[row + 1];
 		while (rowStart < rowEnd) {
-			// If the lement is found, write the content to the specified
-			// position and return true.
+			// If the lement is found, return it.
 			if (columnIndications[rowStart] == col) {
 				return valueStorage[rowStart];
 			}
+			
 			// If the column of the current element is already larger than the
 			// requested column, the requested element cannot be contained
 			// in the matrix and we may therefore stop searching.
@@ -370,6 +417,8 @@ public:
 
 	/*!
 	 * Returns the number of rows of the matrix.
+	 *
+	 * @returns The number of rows of the matrix.
 	 */
 	uint_fast64_t getRowCount() const {
 		return rowCount;
@@ -377,43 +426,18 @@ public:
 
 	/*!
 	 * Returns the number of columns of the matrix.
+	 *
+	 * @returns The number of columns of the matrix.
 	 */
 	uint_fast64_t getColumnCount() const {
 		return colCount;
 	}
 
-	/*!
-	 * Returns a pointer to the value storage of the matrix.
-	 * @return A pointer to the value storage of the matrix.
-	 */
-	std::vector<T> const& getStorage() const {
-		return valueStorage;
-	}
-
-	/*!
-	 * Returns a pointer to the array that stores the start indices of non-zero
-	 * entries in the value storage for each row.
-	 * @return A pointer to the array that stores the start indices of non-zero
-	 * entries in the value storage for each row.
-	 */
-	std::vector<uint_fast64_t> const& getRowIndications() const {
-		return rowIndications;
-	}
-
-	/*!
-	 * Returns a pointer to an array that stores the column of each non-zero
-	 * element.
-	 * @return A pointer to an array that stores the column of each non-zero
-	 * element.
-	 */
-	std::vector<uint_fast64_t> const& getColumnIndications() const {
-		return columnIndications;
-	}
-
 	/*!
 	 * Checks whether the internal status of the matrix makes it ready for
 	 * reading access.
-	 * @return True iff the internal status of the matrix makes it ready for
+	 *
+	 * @returns True iff the internal status of the matrix makes it ready for
 	 * reading access.
 	 */
 	bool isReadReady() {
@@ -423,7 +447,8 @@ public:
 	/*!
 	 * Checks whether the matrix was initialized previously. The matrix may
 	 * still require to be finalized, even if this check returns true.
-	 * @return True iff the matrix was initialized previously.
+	 *
+	 * @returns True iff the matrix was initialized previously.
 	 */
 	bool isInitialized() {
 		return (internalStatus == MatrixStatus::Initialized || internalStatus == MatrixStatus::ReadReady);
@@ -431,7 +456,8 @@ public:
 
 	/*!
 	 * Returns the internal state of the matrix.
-	 * @return The internal state of the matrix.
+	 *
+	 * @returns The internal state of the matrix.
 	 */
 	MatrixStatus getState() {
 		return internalStatus;
@@ -439,15 +465,17 @@ public:
 
 	/*!
 	 * Checks whether the internal state of the matrix signals an error.
-	 * @return True iff the internal state of the matrix signals an error.
+	 *
+	 * @returns True iff the internal state of the matrix signals an error.
 	 */
 	bool hasError() const {
 		return (internalStatus == MatrixStatus::Error);
 	}
 
 	/*!
-	 * Returns the number of non-zero entries that are not on the diagonal.
-	 * @returns The number of non-zero entries that are not on the diagonal.
+	 * Returns the number of non-zero entries in the matrix.
+	 *
+	 * @returns The number of non-zero entries in the matrix.
 	 */
 	uint_fast64_t getNonZeroEntryCount() const {
 		return nonZeroEntryCount;
@@ -455,40 +483,42 @@ public:
 
 	/*!
 	 * This function makes the rows given by the bit vector absorbing.
+	 *
 	 * @param rows A bit vector indicating which rows to make absorbing.
-	 * @return True iff the operation was successful.
+	 * @returns True iff the operation was successful.
 	 */
-	bool makeRowsAbsorbing(const storm::storage::BitVector rows) {
+	bool makeRowsAbsorbing(storm::storage::BitVector const& rows) {
 		bool result = true;
 		for (auto row : rows) {
 			result &= makeRowAbsorbing(row, row);
 		}
-
 		return result;
 	}
 
 	/*!
 	 * This function makes the groups of rows given by the bit vector absorbing.
-	 * @param rows A bit vector indicating which row groups to make absorbing.
+	 *
+	 * @param rowGroupConstraint A bit vector indicating which row groups to make absorbing.
+	 * @param rowGroupIndices A vector indicating which rows belong to a given row group.
 	 * @return True iff the operation was successful.
 	 */
-	bool makeRowsAbsorbing(const storm::storage::BitVector rows, std::vector<uint_fast64_t> const& nondeterministicChoices) {
+	bool makeRowsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices) {
 		bool result = true;
-		for (auto index : rows) {
-			for (uint_fast64_t row = nondeterministicChoices[index]; row < nondeterministicChoices[index + 1]; ++row) {
-				result &= makeRowAbsorbing(row, index);
+		for (auto rowGroup : rowGroupConstraint) {
+			for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) {
+				result &= makeRowAbsorbing(row, rowGroup);
 			}
 		}
-
 		return result;
 	}
 
 
 	/*!
-	 * This function makes the given row absorbing. This means that all
-	 * entries in will be set to 0 and the value 1 will be written
-	 * to the element on the (pseudo-) diagonal.
+	 * This function makes the given row absorbing. This means that all entries will be set to 0
+	 * except the one at the specified column, which is set to 1 instead.
+	 *
 	 * @param row The row to be made absorbing.
+	 * @param column The index of the column whose value is to be set to 1.
 	 * @returns True iff the operation was successful.
 	 */
 	bool makeRowAbsorbing(const uint_fast64_t row, const uint_fast64_t column) {
@@ -504,14 +534,18 @@ public:
 		uint_fast64_t rowStart = rowIndications[row];
 		uint_fast64_t rowEnd = rowIndications[row + 1];
 
+		// If the row has no elements in it, we cannot make it absorbing, because we would need to
+		// move all elements in the vector of nonzeros otherwise.
 		if (rowStart >= rowEnd) {
 			LOG4CPLUS_ERROR(logger, "Cannot make row " << row << " absorbing, because there is no entry in this row.");
 			throw storm::exceptions::InvalidStateException() << "Cannot make row " << row << " absorbing, because there is no entry in this row.";
 		}
 
+		// If there is at least one nonzero entry in this row, we can just set it to one, modify its
+		// column indication to the one given by the parameter and set all subsequent elements of this
+		// row to zero.
 		valueStorage[rowStart] = storm::utility::constGetOne<T>();
 		columnIndications[rowStart] = column;
-
 		for (uint_fast64_t index = rowStart + 1; index < rowEnd; ++index) {
 			valueStorage[index] = storm::utility::constGetZero<T>();
 			columnIndications[index] = 0;
@@ -521,14 +555,15 @@ public:
 	}
 
 	/*
-	 * Computes the sum of the elements in the given row whose column bits
-	 * are set to one on the given constraint.
+	 * Computes the sum of the elements in the given row whose column bits are set to one on the
+	 * given constraint.
+	 *
 	 * @param row The row whose elements to add.
 	 * @param constraint A bit vector that indicates which columns to add.
-	 * @return The sum of the elements in the given row whose column bits
+	 * @returns The sum of the elements in the given row whose column bits
 	 * are set to one on the given constraint.
 	 */
-	T getConstrainedRowSum(const uint_fast64_t row, const storm::storage::BitVector& constraint) const {
+	T getConstrainedRowSum(uint_fast64_t row, storm::storage::BitVector const& constraint) const {
 		T result(0);
 		for (uint_fast64_t i = rowIndications[row]; i < rowIndications[row + 1]; ++i) {
 			if (constraint.get(columnIndications[i])) {
@@ -539,44 +574,55 @@ public:
 	}
 
 	/*!
-	 * Computes a vector in which each element is the sum of those elements in the
-	 * corresponding row whose column bits are set to one in the given constraint.
+	 * Computes a vector whose ith element is the sum of the elements in the ith row where only
+	 * those elements are added whose bits are set to true in the given column constraint and only
+	 * those rows are treated whose bits are set to true in the given row constraint.
+	 *
 	 * @param rowConstraint A bit vector that indicates for which rows to perform summation.
 	 * @param columnConstraint A bit vector that indicates which columns to add.
-	 * @param resultVector A pointer to the resulting vector that has at least
-	 * as many elements as there are bits set to true in the constraint.
+	 *
+	 * @returns A vector whose ith element is the sum of the elements in the ith row where only
+	 * those elements are added whose bits are set to true in the given column constraint and only
+	 * those rows are treated whose bits are set to true in the given row constraint.
 	 */
-	void getConstrainedRowSumVector(const storm::storage::BitVector& rowConstraint, const storm::storage::BitVector& columnConstraint, std::vector<T>* resultVector) const {
+	std::vector<T> getConstrainedRowSumVector(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint) const {
+		std::vector<T> result(rowConstraint.getNumberOfSetBits());
 		uint_fast64_t currentRowCount = 0;
 		for (auto row : rowConstraint) {
-			(*resultVector)[currentRowCount++] = getConstrainedRowSum(row, columnConstraint);
+			result[currentRowCount++] = getConstrainedRowSum(row, columnConstraint);
 		}
+		return result;
 	}
 
 	/*!
-	 * Computes a vector in which each element is the sum of those elements in the
-	 * corresponding row whose column bits are set to one in the given constraint.
-	 * @param rowConstraint A bit vector that indicates for which rows to perform summation.
+	 * Computes a vector whose elements represent the sums of selected (given by the column
+	 * constraint) entries for all rows in selected row groups given by the row group constraint.
+	 *
+	 * @param rowGroupConstraint A bit vector that indicates which row groups are to be considered.
+	 * @param rowGroupIndices A vector indicating which rows belong to a given row group.
 	 * @param columnConstraint A bit vector that indicates which columns to add.
-	 * @param resultVector A pointer to the resulting vector that has at least
-	 * as many elements as there are bits set to true in the constraint.
+	 * @returns 
 	 */
-	void getConstrainedRowSumVector(const storm::storage::BitVector& rowConstraint, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, const storm::storage::BitVector& columnConstraint, std::vector<T>* resultVector) const {
+	std::vector<T> getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, storm::storage::BitVector const& columnConstraint, uint_fast64_t numberOfRows) const {
+		std::vector<T> result(numberOfRows);
 		uint_fast64_t currentRowCount = 0;
-		for (auto index : rowConstraint) {
-			for (uint_fast64_t row = nondeterministicChoiceIndices[index]; row < nondeterministicChoiceIndices[index + 1]; ++row) {
-				(*resultVector)[currentRowCount++] = getConstrainedRowSum(row, columnConstraint);
+		for (auto rowGroup : rowGroupConstraint) {
+			for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) {
+				result[currentRowCount++] = getConstrainedRowSum(row, columnConstraint);
 			}
 		}
+		return result;
 	}
 
 	/*!
-	 * Creates a sub-matrix of the current matrix by dropping all rows and
-	 * columns whose bits are not set to one in the given bit vector.
-	 * @param constraint A bit vector indicating which rows and columns to drop.
-	 * @return A pointer to a sparse matrix that is a sub-matrix of the current one.
+	 * Creates a submatrix of the current matrix by dropping all rows and columns whose bits are not
+	 * set to one in the given bit vector.
+	 *
+	 * @param constraint A bit vector indicating which rows and columns to keep.
+	 * @returns A matrix corresponding to a submatrix of the current matrix in which only rows and
+	 * columns given by the constraint are kept and all others are dropped.
 	 */
-	SparseMatrix* getSubmatrix(storm::storage::BitVector& constraint) const {
+	SparseMatrix getSubmatrix(storm::storage::BitVector const& constraint) const {
 		LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix with " << constraint.getNumberOfSetBits() << " rows.");
 
 		// Check for valid constraint.
@@ -597,8 +643,8 @@ public:
 		}
 
 		// Create and initialize resulting matrix.
-		SparseMatrix* result = new SparseMatrix(constraint.getNumberOfSetBits());
-		result->initialize(subNonZeroEntries);
+		SparseMatrix result(constraint.getNumberOfSetBits());
+		result.initialize(subNonZeroEntries);
 
 		// Create a temporary array that stores for each index whose bit is set
 		// to true the number of bits that were set before that particular index.
@@ -617,7 +663,7 @@ public:
 		for (auto rowIndex : constraint) {
 			for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) {
 				if (constraint.get(columnIndications[i])) {
-					result->addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[i]], valueStorage[i]);
+					result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[i]], valueStorage[i]);
 				}
 			}
 
@@ -628,47 +674,50 @@ public:
 		delete[] bitsSetBeforeIndex;
 
 		// Finalize sub-matrix and return result.
-		result->finalize();
+		result.finalize();
 		LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix.");
 		return result;
 	}
 
-	SparseMatrix* getSubmatrix(storm::storage::BitVector& constraint, std::vector<uint_fast64_t> const& rowIndices) const {
+	/*!
+	 * Creates a submatrix of the current matrix by keeping only row groups and columns in the given
+	 * row group constraint.
+	 *
+	 * @param rowGroupConstraint A bit vector indicating which row groups and columns to keep.
+	 * @param rowGroupIndices A vector indicating which rows belong to a given row group.
+	 * @returns A matrix corresponding to a submatrix of the current matrix in which only row groups
+	 * and columns given by the row group constraint are kept and all others are dropped.
+	 */
+	SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices) const {
 		LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size).");
 
-		// Check for valid constraint.
-		if (constraint.getNumberOfSetBits() == 0) {
-			LOG4CPLUS_ERROR(logger, "Trying to create a sub-matrix of size 0.");
-			throw storm::exceptions::InvalidArgumentException("Trying to create a sub-matrix of size 0.");
-		}
-
 		// First, we need to determine the number of non-zero entries and the number of rows of the
 		// sub-matrix.
 		uint_fast64_t subNonZeroEntries = 0;
 		uint_fast64_t subRowCount = 0;
-		for (auto index : constraint) {
-			subRowCount += rowIndices[index + 1] - rowIndices[index];
-			for (uint_fast64_t i = rowIndices[index]; i < rowIndices[index + 1]; ++i) {
+		for (auto index : rowGroupConstraint) {
+			subRowCount += rowGroupIndices[index + 1] - rowGroupIndices[index];
+			for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) {
 				for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) {
-					if (constraint.get(columnIndications[j])) {
+					if (rowGroupConstraint.get(columnIndications[j])) {
 						++subNonZeroEntries;
 					}
 				}
 			}
 		}
 
-		LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << constraint.getNumberOfSetBits() << ".");
+		LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << rowGroupConstraint.getNumberOfSetBits() << ".");
 
 		// Create and initialize resulting matrix.
-		SparseMatrix* result = new SparseMatrix(subRowCount, constraint.getNumberOfSetBits());
-		result->initialize(subNonZeroEntries);
+		SparseMatrix result(subRowCount, rowGroupConstraint.getNumberOfSetBits());
+		result.initialize(subNonZeroEntries);
 
 		// Create a temporary array that stores for each index whose bit is set
 		// to true the number of bits that were set before that particular index.
 		uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[colCount];
 		uint_fast64_t lastIndex = 0;
 		uint_fast64_t currentNumberOfSetBits = 0;
-		for (auto index : constraint) {
+		for (auto index : rowGroupConstraint) {
 			while (lastIndex <= index) {
 				bitsSetBeforeIndex[lastIndex++] = currentNumberOfSetBits;
 			}
@@ -677,11 +726,11 @@ public:
 
 		// Copy over selected entries.
 		uint_fast64_t rowCount = 0;
-		for (auto index : constraint) {
-			for (uint_fast64_t i = rowIndices[index]; i < rowIndices[index + 1]; ++i) {
+		for (auto index : rowGroupConstraint) {
+			for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) {
 				for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) {
-					if (constraint.get(columnIndications[j])) {
-						result->addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]);
+					if (rowGroupConstraint.get(columnIndications[j])) {
+						result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]);
 					}
 				}
 				++rowCount;
@@ -692,11 +741,16 @@ public:
 		delete[] bitsSetBeforeIndex;
 
 		// Finalize sub-matrix and return result.
-		result->finalize();
+		result.finalize();
 		LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix.");
 		return result;
 	}
 
+	/*!
+	 * Performs a change to the matrix that is needed if this matrix is to be interpreted as a
+	 * set of linear equations. In particular, it transforms A to (1-A), meaning that the elements
+	 * on the diagonal are inverted with respect to addition and the other elements are negated.
+	 */
 	void convertToEquationSystem() {
 		invertDiagonal();
 		negateAllNonDiagonalElements();
@@ -704,28 +758,34 @@ public:
 
 	/*!
 	 * Inverts all elements on the diagonal, i.e. sets the diagonal values to 1 minus their previous
-	 * value.
-	 * Requires the matrix to contain each diagonal element AND to be square!
+	 * value. Requires the matrix to contain each diagonal element and to be square.
 	 */
 	void invertDiagonal() {
+		// Check if the matrix is square, because only then it makes sense to perform this
+		// transformation.
 		if (this->getRowCount() != this->getColumnCount()) {
 			throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to be square!";
 		}
+		
+		// Now iterate over all rows and set the diagonal elements to the inverted value.
+		// If there is a row without the diagonal element, an exception is thrown.
 		T one = storm::utility::constGetOne<T>();
-		bool foundRow;
+		bool foundDiagonalElement = false;
 		for (uint_fast64_t row = 0; row < rowCount; ++row) {
 			uint_fast64_t rowStart = rowIndications[row];
 			uint_fast64_t rowEnd = rowIndications[row + 1];
-			foundRow = false;
+			foundDiagonalElement = false;
 			while (rowStart < rowEnd) {
 				if (columnIndications[rowStart] == row) {
 					valueStorage[rowStart] = one - valueStorage[rowStart];
-					foundRow = true;
+					foundDiagonalElement = true;
 					break;
 				}
 				++rowStart;
 			}
-			if (!foundRow) {
+			
+			// Throw an exception if a row did not have an element on the diagonal.
+			if (!foundDiagonalElement) {
 				throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to contain all diagonal entries!";
 			}
 		}
@@ -735,9 +795,13 @@ public:
 	 * Negates all non-zero elements that are not on the diagonal.
 	 */
 	void negateAllNonDiagonalElements() {
+		// Check if the matrix is square, because only then it makes sense to perform this
+		// transformation.
 		if (this->getRowCount() != this->getColumnCount()) {
 			throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to be square!";
 		}
+		
+		// Iterate over all rows and negate all the elements that are not on the diagonal.
 		for (uint_fast64_t row = 0; row < rowCount; ++row) {
 			uint_fast64_t rowStart = rowIndications[row];
 			uint_fast64_t rowEnd = rowIndications[row + 1];
@@ -759,7 +823,7 @@ public:
 		uint_fast64_t rowCount = this->getRowCount();
 		uint_fast64_t colCount = this->getColumnCount();
 		if (rowCount != colCount) {
-			throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::getJacobiDecomposition requires the Matrix to be square!";
+			throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::getJacobiDecomposition requires the Matrix to be square.";
 		}
 		storm::storage::SparseMatrix<T> *resultLU = new storm::storage::SparseMatrix<T>(*this);
 		storm::storage::SparseMatrix<T> *resultDinv = new storm::storage::SparseMatrix<T>(rowCount, colCount);
@@ -781,15 +845,16 @@ public:
 	/*!
 	 * Performs a pointwise matrix multiplication of the matrix with the given matrix and returns a
 	 * vector containing the sum of the elements in each row of the resulting matrix.
+	 *
 	 * @param otherMatrix A reference to the matrix with which to perform the pointwise multiplication.
 	 * This matrix must be a submatrix of the current matrix in the sense that it may not have
 	 * non-zero entries at indices where there is a zero in the current matrix.
-	 * @return A vector containing the sum of the elements in each row of the matrix resulting from
+	 * @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) {
 		// Prepare result.
-		std::vector<T>* result = new std::vector<T>(rowCount, storm::utility::constGetZero<T>());
+		std::vector<T> result(rowCount, storm::utility::constGetZero<T>());
 
 		// Iterate over all elements of the current matrix and either continue with the next element
 		// in case the given matrix does not have a non-zero element at this column position, or
@@ -802,7 +867,7 @@ public:
 					// If the precondition of this method (i.e. that the given matrix is a submatrix
 					// of the current one) was fulfilled, we know now that the two elements are in
 					// the same column, so we can multiply and add them to the row sum vector.
-					(*result)[row] += otherMatrix.valueStorage[nextOtherElement] * valueStorage[element];
+					result[row] += otherMatrix.valueStorage[nextOtherElement] * valueStorage[element];
 					++nextOtherElement;
 				}
 			}
@@ -810,84 +875,39 @@ public:
 
 		return result;
 	}
-
-	T multiplyRowWithVector(constRowsIterator& rowsIt, uint_fast64_t rowsIte, std::vector<T>& vector) const {
-		T result = storm::utility::constGetZero<T>();
-		for (; rowsIt != rowsIte; ++rowsIt) {
-			result += (rowsIt.value()) * vector[rowsIt.column()];
-		}
-		return result;
-	}
-
-	void multiplyWithVector(std::vector<T>& vector, std::vector<T>& result) const {
-		typename std::vector<T>::iterator resultIt = result.begin();
-		typename std::vector<T>::iterator resultIte = result.end();
-		constRowsIterator rowIt = this->constRowsIteratorBegin();
-		uint_fast64_t nextRow = 1;
-
-		for (; resultIt != resultIte; ++resultIt, ++nextRow) {
-			*resultIt = multiplyRowWithVector(rowIt, this->rowIndications[nextRow], vector);
-		}
-	}
-
-	void multiplyWithVectorInPlace(std::vector<T>& vector) const {
-		typename std::vector<T>::iterator resultIt = vector.begin();
-		typename std::vector<T>::iterator resultIte = vector.end();
-		constRowsIterator rowIt = this->constRowsIteratorBegin();
-		uint_fast64_t nextRow = 1;
-
-		for (; resultIt != resultIte; ++resultIt, ++nextRow) {
-			*resultIt = multiplyRowWithVector(rowIt, this->rowIndications[nextRow], vector);
-		}
-	}
-
-	void multiplyWithVector(std::vector<uint_fast64_t> const& states, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<T>& vector, std::vector<T>& result) const {
-		constRowsIterator rowsIt = this->constRowsIteratorBegin();
-		uint_fast64_t nextRow = 1;
-
-		for (auto stateIt = states.cbegin(), stateIte = states.cend(); stateIt != stateIte; ++stateIt) {
-			rowsIt.setOffset(this->rowIndications[nondeterministicChoiceIndices[*stateIt]]);
-			nextRow = nondeterministicChoiceIndices[*stateIt] + 1;
-			for (auto rowIt = nondeterministicChoiceIndices[*stateIt], rowIte = nondeterministicChoiceIndices[*stateIt + 1]; rowIt != rowIte; ++rowIt, ++nextRow) {
-				result[rowIt] = multiplyRowWithVector(rowsIt, this->rowIndications[nextRow], vector);
-			}
-		}
-	}
-
-	void multiplyWithVectorInPlace(std::vector<uint_fast64_t> const& states, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<T>& vector) const {
-		constRowsIterator rowsIt = this->constRowsIteratorBegin();
-		uint_fast64_t nextRow = 1;
-
-		for (auto stateIt = states.cbegin(), stateIte = states.cend(); stateIt != stateIte; ++stateIt) {
-			rowsIt.setOffset(this->rowIndications[nondeterministicChoiceIndices[*stateIt]]);
-			nextRow = nondeterministicChoiceIndices[*stateIt] + 1;
-			for (auto rowIt = nondeterministicChoiceIndices[*stateIt], rowIte = nondeterministicChoiceIndices[*stateIt + 1]; rowIt != rowIte; ++rowIt, ++nextRow) {
-				vector[rowIt] = multiplyRowWithVector(rowsIt, this->rowIndications[nextRow], vector);
-			}
-		}
-	}
-
-	void multiplyAddAndReduceInPlace(std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<T>& vector, std::vector<T> const& summand, bool minimize) const {
-		constRowsIterator rowsIt = this->constRowsIteratorBegin();
-		uint_fast64_t nextRow = 1;
-
-		for (uint_fast64_t stateIt = 0, stateIte = nondeterministicChoiceIndices.size() - 1; stateIt < stateIte; ++stateIt) {
-			vector[stateIt] = multiplyRowWithVector(rowsIt, this->rowIndications[nextRow], vector) + summand[nondeterministicChoiceIndices[stateIt]];
-			++nextRow;
-			for (uint_fast64_t rowIt = nondeterministicChoiceIndices[stateIt] + 1, rowIte = nondeterministicChoiceIndices[stateIt + 1]; rowIt != rowIte; ++rowIt, ++nextRow) {
-				T value = multiplyRowWithVector(rowsIt, this->rowIndications[nextRow], vector) + summand[rowIt];
-				if (minimize && value < vector[stateIt]) {
-					vector[stateIt] = value;
-				} else if (!minimize && value > vector[stateIt]) {
-					vector[stateIt] = value;
-				}
+	
+	/*!
+	 * Multiplies the matrix with the given vector and writes the result to given result vector.
+	 *
+	 * @param vector The vector with which to multiply the matrix.
+	 * @param result The vector that is supposed to hold the result of the multiplication after the
+	 * operation.
+	 * @returns The product of the matrix and the given vector as the content of the given result
+	 * vector.
+	 */
+	void multiplyWithVector(std::vector<T> const& vector, std::vector<T>& result) const {
+		// Initialize two iterators that 
+		ConstRowsIterator elementIt(*this);
+		ConstRowsIterator elementIte(*this, 1);
+		
+		// Iterate over all positions of the result vector and compute its value as the scalar
+		// product of the corresponding row with the input vector.
+		// Note that only the end iterator has to be moved by one row, because the other iterator
+		// is automatically moved forward one row by the inner loop.
+		for (uint_fast64_t i = 0; i < result.size(); ++i, elementIte.moveToNextRow()) {
+			result[i] = storm::utility::constGetZero<T>();
+			
+			// Perform the scalar product.
+			for (; elementIt != elementIte; ++elementIt) {
+				result[i] += elementIt.value() * vector[elementIt.column()];
 			}
 		}
 	}
 
 	/*!
 	 * Returns the size of the matrix in memory measured in bytes.
-	 * @return The size of the matrix in memory measured in bytes.
+	 *
+	 * @returns The size of the matrix in memory measured in bytes.
 	 */
 	uint_fast64_t getSizeInMemory() const {
 		uint_fast64_t size = sizeof(*this);
@@ -900,87 +920,122 @@ public:
 		return size;
 	}
 
-	constRowsIterator constRowsIteratorBegin() const {
-		return constRowsIterator(*this);
+	/*!
+	 * Returns a const iterator to the elements of the matrix.
+	 * @param row If given, this specifies the starting row of the iterator.
+	 *
+	 * @returns An iterator to the elements of the matrix that cannot be used for modifying the
+	 * contents.
+	 */
+	ConstRowsIterator begin(uint_fast64_t row = 0) const {
+		return ConstRowsIterator(*this, row);
+	}
+	
+	/*!
+	 * Returns a const iterator that points to the first element after the matrix.
+	 *
+	 * @returns A const iterator that points past the last element of the matrix.
+	 */
+	ConstRowsIterator end() const {
+		return ConstRowsIterator(this->rowCount);
+	}
+	
+	/*!
+	 * Returns a const iterator that points to the first element after the given row.
+	 *
+	 * @returns row The row past which this iterator has to point.
+	 * @returns A const iterator that points to the first element after the given row.
+	 */
+	ConstRowsIterator end(uint_fast64_t row) const {
+		return ConstRowsIterator(row + 1);
 	}
 
 	/*!
-	 * Returns an iterator to the columns of the non-zero entries of the given
-	 * row.
-	 * @param row The row whose columns the iterator will return.
-	 * @return An iterator to the columns of the non-zero entries of the given
-	 * row.
+	 * Returns an iterator to the columns of the non-zero entries of the matrix.
+	 *
+	 * @param row If given, the iterator will start at the specified row.
+	 * @returns An iterator to the columns of the non-zero entries of the matrix.
 	 */
-	constIndexIterator beginConstColumnIterator(uint_fast64_t row) const {
+	ConstIndexIterator constColumnIteratorBegin(uint_fast64_t row = 0) const {
 		return &(this->columnIndications[0]) + this->rowIndications[row];
 	}
 
 	/*!
-	 * Returns an iterator referring to the element after the given row.
-	 * @param row The row for which the iterator should point to the past-the-end
-	 * element.
+	 * Returns an iterator that points to the first element after the matrix.
+	 *
+	 * @returns An iterator that points to the first element after the matrix.
 	 */
-	constIndexIterator endConstColumnIterator(uint_fast64_t row) const {
+	ConstIndexIterator constColumnIteratorEnd() const {
+		return &(this->columnIndications[0]) + this->rowIndications[rowCount];
+	}
+	
+	/*!
+	 * Returns an iterator that points to the first element after the given row.
+	 *
+	 * @param row The row past which this iterator has to point.
+	 * @returns An iterator that points to the first element after the matrix.
+	 */
+	ConstIndexIterator constColumnIteratorEnd(uint_fast64_t row) const {
 		return &(this->columnIndications[0]) + this->rowIndications[row + 1];
 	}
-
+	
 	/*!
-	 *	Returns an iterator over the elements of the given row. The iterator
-	 *	will include no zero entries.
-	 *	@param row The row whose elements the iterator will return.
-	 *	@return An iterator over the elements of the given row.
+	 * Returns an iterator to the values of the non-zero entries of the matrix.
+	 *
+	 * @param row If given, the iterator will start at the specified row.
+	 * @returns An iterator to the values of the non-zero entries of the matrix.
 	 */
-	constIterator beginConstIterator(uint_fast64_t row) const {
+	ConstValueIterator constValueIteratorBegin(uint_fast64_t row = 0) const {
 		return &(this->valueStorage[0]) + this->rowIndications[row];
 	}
+	
 	/*!
-	 *	Returns an iterator pointing to the first element after the given
-	 *	row.
-	 *	@param row The row for which the iterator should point to the
-	 *	past-the-end element.
-	 *	@return An iterator to the element after the given row.
+	 * Returns an iterator that points to the first element after the matrix.
+	 *
+	 * @returns An iterator that points to the first element after the matrix.
 	 */
-	constIterator endConstIterator(uint_fast64_t row) const {
+	ConstValueIterator constValueIteratorEnd() const {
+		return &(this->valueStorage[0]) + this->rowIndications[rowCount];
+	}
+	
+	/*!
+	 * Returns an iterator that points to the first element after the given row.
+	 *
+	 * @param row The row past which this iterator has to point.
+	 * @returns An iterator that points to the first element after the matrix.
+	 */
+	ConstValueIterator constValueIteratorEnd(uint_fast64_t row) const {
 		return &(this->valueStorage[0]) + this->rowIndications[row + 1];
 	}
-
+	
 	/*!
-	 *	@brief Calculate sum of all entries in given row.
+	 * Computes the sum of the elements in a given row.
 	 *
-	 *	Adds up all values in the given row
-	 *	and returns the sum.
-	 *	@param row The row that should be added up.
-	 *	@return Sum of the row.
+	 * @param row The row that should be summed.
+	 * @return Sum of the row.
 	 */
 	T getRowSum(uint_fast64_t row) const {
 		T sum = storm::utility::constGetZero<T>();
-		for (auto it = this->beginConstIterator(row); it != this->endConstIterator(row); it++) {
+		for (auto it = this->constValueIteratorBegin(row), ite = this->constValueIteratorEnd(row); it != ite; ++it) {
 			sum += *it;
 		}
 		return sum;
 	}
 
 	/*!
-	 *	@brief Checks if it is a submatrix of the given matrix.
+	 * Checks if the given matrix is a submatrix of the current matrix, where A matrix A is called a
+	 * submatrix of B if a value in A is only nonzero, if the value in B at the same position is
+	 * also nonzero. Furthermore, A and B have to have the same size.
 	 *
-	 *	A matrix A is a submatrix of B if a value in A is only nonzero, if
-	 *	the value in B at the same position is also nonzero. Furthermore, A
-	 *	and B have to have the same size.
-	 *	@param matrix Matrix to check against.
-	 *	@return True iff this is a submatrix of matrix.
+	 * @param matrix The matrix that is possibly a submatrix of the current matrix.
+	 * @returns True iff the given matrix is a submatrix of the current matrix.
 	 */
-	bool isSubmatrixOf(SparseMatrix<T> const& matrix) const {
+	bool containsAllPositionsOf(SparseMatrix<T> const& matrix) const {
+		// Check for mismatching sizes.
 		if (this->getRowCount() != matrix.getRowCount()) return false;
 		if (this->getColumnCount() != matrix.getColumnCount()) return false;
 
-		/*
-		for (uint_fast64_t row = 0; row < this->getRowCount(); ++row) {
-			for (uint_fast64_t elem = rowIndications[row], elem2 = matrix.rowIndications[row]; elem < rowIndications[row + 1] && elem < matrix.rowIndications[row + 1]; ++elem, ++elem2) {
-				if (columnIndications[elem] < matrix.columnIndications[elem2]) return false;
-			}
-		}
-		*/
-
+		// Check the subset property for all rows individually.
 		for (uint_fast64_t row = 0; row < this->getRowCount(); ++row) {
 			for (uint_fast64_t elem = rowIndications[row], elem2 = matrix.rowIndications[row]; elem < rowIndications[row + 1] && elem < matrix.rowIndications[row + 1]; ++elem) {
 				// Skip over all entries of the other matrix that are before the current entry in
@@ -996,7 +1051,8 @@ public:
 
 	/*!
 	 * Retrieves a compressed string representation of the matrix.
-	 * @return a compressed string representation of the matrix.
+	 *
+	 * @returns a compressed string representation of the matrix.
 	 */
 	std::string toStringCompressed() const {
 		std::stringstream result;
@@ -1008,13 +1064,14 @@ public:
 
 	/*!
 	 * Retrieves a (non-compressed) string representation of the matrix.
-	 * Note: the matrix is presented densely. That is, all zeros are filled in and are part of the string
-	 * representation.
-	 * @param nondeterministicChoiceIndices A vector indicating which rows belong together. If given, rows belonging
-	 * to separate groups will be separated by a dashed line.
-	 * @return a (non-compressed) string representation of the matrix.
+	 * Note: the matrix is presented densely. That is, all zeros are filled in and are part of the
+	 * string representation, so calling this method on big matrices should be done with care.
+	 *
+	 * @param rowGroupIndices A vector indicating to which group of rows a given row belongs. If
+	 * given, rows of different groups will be separated by a dashed line.
+	 * @returns A (non-compressed) string representation of the matrix.
 	 */
-	std::string toString(std::vector<uint_fast64_t> const* nondeterministicChoiceIndices) const {
+	std::string toString(std::vector<uint_fast64_t> const* rowGroupIndices = nullptr) const {
 		std::stringstream result;
 		uint_fast64_t currentNondeterministicChoiceIndex = 0;
 
@@ -1029,9 +1086,9 @@ public:
 		for (uint_fast64_t i = 0; i < rowCount; ++i) {
 			uint_fast64_t nextIndex = rowIndications[i];
 
-			// If we need to group of rows, print a dashed line in case we have moved to the next group of rows.
-			if (nondeterministicChoiceIndices != nullptr) {
-				if (i == (*nondeterministicChoiceIndices)[currentNondeterministicChoiceIndex]) {
+			// If we need to group rows, print a dashed line in case we have moved to the next group of rows.
+			if (rowGroupIndices != nullptr) {
+				if (i == (*rowGroupIndices)[currentNondeterministicChoiceIndex]) {
 					if (i != 0) {
 						result << "\t(\t";
 						for (uint_fast64_t j = 0; j < colCount - 2; ++j) {
@@ -1140,45 +1197,36 @@ private:
 	}
 
 	/*!
-	 * Prepares the internal CSR storage. For this, it requires
-	 * non_zero_entry_count and row_count to be set correctly.
-	 * @param alsoPerformAllocation If set to true, all entries are pre-allocated. This is the default.
+	 * Prepares the internal storage. For this, it requires the number of non-zero entries and the
+	 * amount of rows to be set correctly.
+	 *
+	 * @param initializeElements If set to true, all entries are initialized.
 	 * @return True on success, false otherwise (allocation failed).
 	 */
-	bool prepareInternalStorage(const bool alsoPerformAllocation) {
-		if (alsoPerformAllocation) {
+	bool prepareInternalStorage(bool initializeElements = true) {
+		if (initializeElements) {
 			// Set up the arrays for the elements that are not on the diagonal.
 			valueStorage.resize(nonZeroEntryCount, storm::utility::constGetZero<T>());
 			columnIndications.resize(nonZeroEntryCount, 0);
 
-			// Set up the row_indications vector and reserve one element more than
-			// there are rows in order to put a sentinel element at the end,
-			// which eases iteration process.
+			// Set up the rowIndications vector and reserve one element more than there are rows in
+			// order to put a sentinel element at the end, which eases iteration process.
 			rowIndications.resize(rowCount + 1, 0);
 
 			// Return whether all the allocations could be made without error.
 			return ((valueStorage.capacity() >= nonZeroEntryCount) && (columnIndications.capacity() >= nonZeroEntryCount)
 				&& (rowIndications.capacity() >= (rowCount + 1)));
 		} else {
+			// If it was not requested to initialize the elements, we simply reserve the space.
 			valueStorage.reserve(nonZeroEntryCount);
 			columnIndications.reserve(nonZeroEntryCount);
 			rowIndications.reserve(rowCount + 1);
 			return true;
 		}
 	}
-
-	/*!
-	 * Shorthand for prepareInternalStorage(true)
-	 * @see prepareInternalStorage(const bool)
-	 */
-	bool prepareInternalStorage() {
-		return this->prepareInternalStorage(true);
-	}
-
 };
 
 } // namespace storage
-
 } // namespace storm
 
 #endif // STORM_STORAGE_SPARSEMATRIX_H_
diff --git a/src/storm.cpp b/src/storm.cpp
index 28a95c01c..9f939d6ee 100644
--- a/src/storm.cpp
+++ b/src/storm.cpp
@@ -96,7 +96,7 @@ void printFooter() {
 bool parseOptions(const int argc, const char* argv[]) {
 	storm::settings::Settings* s = nullptr;
 	try {
-		storm::settings::Settings::registerModule<storm::modelChecker::GmmxxDtmcPrctlModelChecker<double>>();
+		storm::settings::Settings::registerModule<storm::modelchecker::GmmxxDtmcPrctlModelChecker<double>>();
 		s = storm::settings::newInstance(argc, argv, nullptr);
 	} catch (storm::exceptions::InvalidSettingsException& e) {
 		std::cout << "Could not recover from settings error: " << e.what() << "." << std::endl;
@@ -140,7 +140,7 @@ void cleanUp() {
 }
 
 void testCheckingDie(storm::models::Dtmc<double>& dtmc) {
-	storm::modelChecker::GmmxxDtmcPrctlModelChecker<double>* mc = new storm::modelChecker::GmmxxDtmcPrctlModelChecker<double>(dtmc);
+	storm::modelchecker::GmmxxDtmcPrctlModelChecker<double>* mc = new storm::modelchecker::GmmxxDtmcPrctlModelChecker<double>(dtmc);
 
 	storm::formula::Ap<double>* oneFormula = new storm::formula::Ap<double>("one");
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(oneFormula);
@@ -177,7 +177,7 @@ void testCheckingCrowds(storm::models::Dtmc<double>& dtmc) {
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(observe0Greater1Formula);
 	storm::formula::ProbabilisticNoBoundOperator<double>* probFormula = new storm::formula::ProbabilisticNoBoundOperator<double>(eventuallyFormula);
 
-	storm::modelChecker::GmmxxDtmcPrctlModelChecker<double>* mc = new storm::modelChecker::GmmxxDtmcPrctlModelChecker<double>(dtmc);
+	storm::modelchecker::GmmxxDtmcPrctlModelChecker<double>* mc = new storm::modelchecker::GmmxxDtmcPrctlModelChecker<double>(dtmc);
 	mc->check(*probFormula);
 	delete probFormula;
 
@@ -203,7 +203,7 @@ void testCheckingSynchronousLeader(storm::models::Dtmc<double>& dtmc, uint_fast6
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(electedFormula);
 	storm::formula::ProbabilisticNoBoundOperator<double>* probFormula = new storm::formula::ProbabilisticNoBoundOperator<double>(eventuallyFormula);
 
-	storm::modelChecker::GmmxxDtmcPrctlModelChecker<double>* mc = new storm::modelChecker::GmmxxDtmcPrctlModelChecker<double>(dtmc);
+	storm::modelchecker::GmmxxDtmcPrctlModelChecker<double>* mc = new storm::modelchecker::GmmxxDtmcPrctlModelChecker<double>(dtmc);
 	mc->check(*probFormula);
 	delete probFormula;
 
@@ -229,7 +229,7 @@ void testCheckingDice(storm::models::Mdp<double>& mdp) {
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(twoFormula);
 	storm::formula::ProbabilisticNoBoundOperator<double>* probFormula = new storm::formula::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
 
-	storm::modelChecker::GmmxxMdpPrctlModelChecker<double>* mc = new storm::modelChecker::GmmxxMdpPrctlModelChecker<double>(mdp);
+	storm::modelchecker::GmmxxMdpPrctlModelChecker<double>* mc = new storm::modelchecker::GmmxxMdpPrctlModelChecker<double>(mdp);
 
 	mc->check(*probFormula);
 	delete probFormula;
@@ -320,7 +320,7 @@ void testCheckingAsynchLeader(storm::models::Mdp<double>& mdp) {
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(electedFormula);
 	storm::formula::ProbabilisticNoBoundOperator<double>* probMinFormula = new storm::formula::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
 
-	storm::modelChecker::GmmxxMdpPrctlModelChecker<double>* mc = new storm::modelChecker::GmmxxMdpPrctlModelChecker<double>(mdp);
+	storm::modelchecker::GmmxxMdpPrctlModelChecker<double>* mc = new storm::modelchecker::GmmxxMdpPrctlModelChecker<double>(mdp);
 
 	mc->check(*probMinFormula);
 	delete probMinFormula;
@@ -368,7 +368,7 @@ void testCheckingConsensus(storm::models::Mdp<double>& mdp) {
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(finishedFormula);
 	storm::formula::ProbabilisticNoBoundOperator<double>* probFormula = new storm::formula::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
 
-	storm::modelChecker::GmmxxMdpPrctlModelChecker<double>* mc = new storm::modelChecker::GmmxxMdpPrctlModelChecker<double>(mdp);
+	storm::modelchecker::GmmxxMdpPrctlModelChecker<double>* mc = new storm::modelchecker::GmmxxMdpPrctlModelChecker<double>(mdp);
 
 	mc->check(*probFormula);
 	delete probFormula;
@@ -450,7 +450,7 @@ void testChecking() {
 		std::shared_ptr<storm::models::Mdp<double>> mdp = parser.getModel<storm::models::Mdp<double>>();
 		mdp->printModelInformationToStream(std::cout);
 
-		// testCheckingDice(*mdp);
+        // testCheckingDice(*mdp);
 		// testCheckingAsynchLeader(*mdp);
 		// testCheckingConsensus(*mdp);
 	} else {
diff --git a/src/utility/GraphAnalyzer.h b/src/utility/GraphAnalyzer.h
index 365e95294..5bbb4295c 100644
--- a/src/utility/GraphAnalyzer.h
+++ b/src/utility/GraphAnalyzer.h
@@ -35,7 +35,7 @@ public:
 	 * probability 1 after the invocation of the function.
 	 */
 	template <typename T>
-	static void performProb01(storm::models::AbstractDeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0, storm::storage::BitVector* statesWithProbability1) {
+	static void performProb01(storm::models::AbstractDeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0, storm::storage::BitVector* statesWithProbability1) {
 		// Check for valid parameters.
 		if (statesWithProbability0 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability0' must not be null.");
@@ -63,7 +63,7 @@ public:
 	 * a positive probability of satisfying phi until psi.
 	 */
 	template <typename T>
-	static void performProbGreater0(storm::models::AbstractDeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbabilityGreater0) {
+	static void performProbGreater0(storm::models::AbstractDeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbabilityGreater0) {
 		// Check for valid parameter.
 		if (statesWithProbabilityGreater0 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbabilityGreater0' must not be null.");
@@ -72,7 +72,7 @@ public:
 
 		// Get the backwards transition relation from the model to ease the search.
 		storm::models::GraphTransitions<T> backwardTransitions(*model.getTransitionMatrix(), false);
-
+        
 		// Add all psi states as the already satisfy the condition.
 		*statesWithProbabilityGreater0 |= psiStates;
 
@@ -110,7 +110,7 @@ public:
 	 * have paths satisfying phi until psi.
 	 */
 	template <typename T>
-	static void performProb1(storm::models::AbstractDeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector const& statesWithProbabilityGreater0, storm::storage::BitVector* statesWithProbability1) {
+	static void performProb1(storm::models::AbstractDeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector const& statesWithProbabilityGreater0, storm::storage::BitVector* statesWithProbability1) {
 		// Check for valid parameter.
 		if (statesWithProbability1 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability1' must not be null.");
@@ -134,7 +134,7 @@ public:
 	 * have paths satisfying phi until psi.
 	 */
 	template <typename T>
-	static void performProb1(storm::models::AbstractDeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability1) {
+	static void performProb1(storm::models::AbstractDeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability1) {
 		// Check for valid parameter.
 		if (statesWithProbability1 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability1' must not be null.");
@@ -149,7 +149,7 @@ public:
 	}
 
 	template <typename T>
-	static void performProb01Max(storm::models::AbstractNondeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0, storm::storage::BitVector* statesWithProbability1) {
+	static void performProb01Max(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0, storm::storage::BitVector* statesWithProbability1) {
 		// Check for valid parameters.
 		if (statesWithProbability0 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability0' must not be null.");
@@ -166,7 +166,7 @@ public:
 	}
 
 	template <typename T>
-	static void performProb0A(storm::models::AbstractNondeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0) {
+	static void performProb0A(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0) {
 		// Check for valid parameter.
 		if (statesWithProbability0 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability0' must not be null.");
@@ -201,7 +201,7 @@ public:
 	}
 
 	template <typename T>
-	static void performProb1E(storm::models::AbstractNondeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability1) {
+	static void performProb1E(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability1) {
 		// Check for valid parameters.
 		if (statesWithProbability1 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability1' must not be null.");
@@ -236,7 +236,7 @@ public:
 						// nondeterminstic choices.
 						for (uint_fast64_t row = (*nondeterministicChoiceIndices)[*it]; row < (*nondeterministicChoiceIndices)[*it + 1]; ++row) {
 							bool allSuccessorsInCurrentStates = true;
-							for (auto colIt = transitionMatrix->beginConstColumnIterator(row); colIt != transitionMatrix->endConstColumnIterator(row); ++colIt) {
+							for (auto colIt = transitionMatrix->constColumnIteratorBegin(row); colIt != transitionMatrix->constColumnIteratorEnd(row); ++colIt) {
 								if (!currentStates->get(*colIt)) {
 									allSuccessorsInCurrentStates = false;
 									break;
@@ -270,7 +270,7 @@ public:
 	}
 
 	template <typename T>
-	static void performProb01Min(storm::models::AbstractNondeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0, storm::storage::BitVector* statesWithProbability1) {
+	static void performProb01Min(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0, storm::storage::BitVector* statesWithProbability1) {
 		// Check for valid parameters.
 		if (statesWithProbability0 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability0' must not be null.");
@@ -287,7 +287,7 @@ public:
 	}
 
 	template <typename T>
-	static void performProb0E(storm::models::AbstractNondeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0) {
+	static void performProb0E(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability0) {
 		// Check for valid parameter.
 		if (statesWithProbability0 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability0' must not be null.");
@@ -321,7 +321,7 @@ public:
 					bool addToStatesWithProbability0 = true;
 					for (auto rowIt = nondeterministicChoiceIndices->begin() + *it; rowIt != nondeterministicChoiceIndices->begin() + *it + 1; ++rowIt) {
 						bool hasAtLeastOneSuccessorWithProbabilityGreater0 = false;
-						for (auto colIt = transitionMatrix->beginConstColumnIterator(*rowIt); colIt != transitionMatrix->endConstColumnIterator(*rowIt); ++colIt) {
+						for (auto colIt = transitionMatrix->constColumnIteratorBegin(*rowIt); colIt != transitionMatrix->constColumnIteratorEnd(*rowIt); ++colIt) {
 							if (statesWithProbability0->get(*colIt)) {
 								hasAtLeastOneSuccessorWithProbabilityGreater0 = true;
 								break;
@@ -347,7 +347,7 @@ public:
 	}
 
 	template <typename T>
-	static void performProb1A(storm::models::AbstractNondeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability1) {
+	static void performProb1A(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability1) {
 		// Check for valid parameters.
 		if (statesWithProbability1 == nullptr) {
 			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability1' must not be null.");
@@ -382,7 +382,7 @@ public:
 						// nondeterminstic choices.
 						bool allSuccessorsInCurrentStatesForAllChoices = true;
 						for (uint_fast64_t row = (*nondeterministicChoiceIndices)[*it]; row < (*nondeterministicChoiceIndices)[*it + 1]; ++row) {
-							for (auto colIt = transitionMatrix->beginConstColumnIterator(row); colIt != transitionMatrix->endConstColumnIterator(row); ++colIt) {
+							for (auto colIt = transitionMatrix->constColumnIteratorBegin(row); colIt != transitionMatrix->constColumnIteratorEnd(row); ++colIt) {
 								if (!currentStates->get(*colIt)) {
 									allSuccessorsInCurrentStatesForAllChoices = false;
 									goto afterCheckLoop;
diff --git a/src/utility/IoUtility.cpp b/src/utility/IoUtility.cpp
deleted file mode 100644
index 92d682a7b..000000000
--- a/src/utility/IoUtility.cpp
+++ /dev/null
@@ -1,71 +0,0 @@
-/*
-* IoUtility.cpp
-*
-*  Created on: 17.10.2012
-*      Author: Thomas Heinemann
-*/
-
-#include "src/utility/IoUtility.h"
-
-#include <fstream>
-
-#include "src/parser/DeterministicSparseTransitionParser.h"
-#include "src/parser/AtomicPropositionLabelingParser.h"
-
-namespace storm {
-
-	namespace utility {
-
-		void dtmcToDot(storm::models::Dtmc<double> const &dtmc, std::string filename) {
-			std::shared_ptr<storm::storage::SparseMatrix<double>> matrix(dtmc.getTransitionMatrix());
-			std::ofstream file;
-			file.open(filename);
-
-			file << "digraph dtmc {\n";
-
-			//Specify the nodes and their labels
-			for (uint_fast64_t i = 1; i < dtmc.getNumberOfStates(); i++) {
-				file << "\t" << i << "[label=\"" << i << "\\n{";
-				char komma=' ';
-				std::set<std::string> propositions = dtmc.getPropositionsForState(i);
-				for(auto it = propositions.begin();
-					it != propositions.end();
-					it++) {
-						file << komma << *it;
-						komma=',';
-				}
-
-				file << " }\"];\n";
-
-			}
-
-			for (uint_fast64_t row = 0; row < dtmc.getNumberOfStates(); row++ ) {
-
-				//Then, iterate through the row and write each non-diagonal value into the file
-				for (	auto it = matrix->beginConstColumnIterator(row);
-					it !=  matrix->endConstColumnIterator(row);
-					it++) {
-						double value = 0;
-						matrix->getValue(row,*it,&value);
-						file << "\t" << row << " -> " << *it << " [label=" << value << "]\n";
-				}
-			}
-
-			file << "}\n";
-			file.close();
-		}
-
-		//TODO: Should this stay here or be integrated in the new parser structure?
-		/*storm::models::Dtmc<double>* parseDTMC(std::string const &tra_file, std::string const &lab_file) {
-		storm::parser::DeterministicSparseTransitionParser tp(tra_file);
-		uint_fast64_t node_count = tp.getMatrix()->getRowCount();
-
-		storm::parser::AtomicPropositionLabelingParser lp(node_count, lab_file);
-
-		storm::models::Dtmc<double>* result = new storm::models::Dtmc<double>(tp.getMatrix(), lp.getLabeling());
-		return result;
-		}*/
-
-	}
-
-}
diff --git a/src/utility/IoUtility.h b/src/utility/IoUtility.h
deleted file mode 100644
index 1f6d5b602..000000000
--- a/src/utility/IoUtility.h
+++ /dev/null
@@ -1,45 +0,0 @@
-/*
- * IoUtility.h
- *
- *  Created on: 17.10.2012
- *      Author: Thomas Heinemann
- */
-
-#ifndef STORM_UTILITY_IOUTILITY_H_
-#define STORM_UTILITY_IOUTILITY_H_
-
-#include "src/models/Dtmc.h"
-
-namespace storm  {
-
-namespace utility {
-
-/*!
-    Creates a DOT file which provides the graph of the DTMC.
-
-    Currently, only a version for DTMCs using probabilities of type double is provided.
-    Adaptions for other types may be included later.
-
-    @param  dtmc     The DTMC to output
-    @param  filename The Name of the file to write in. If the file already exists,
-                     it will be overwritten.
-
- */
-void dtmcToDot(storm::models::Dtmc<double> const &dtmc, std::string filename);
-
-/*!
-    Parses a transition file and a labeling file and produces a DTMC out of them.
-    Note that the labeling file may have at most as many nodes as the transition file!
-
-    @param tra_file String containing the location of the transition file (....tra)
-    @param lab_file String containing the location of the labeling file (....lab)
-    @returns The DTMC described by the two files.
-
- */
-//storm::models::Dtmc<double>* parseDTMC(std::string const &tra_file, std::string const &lab_file);
-
-} //namespace utility
-
-} //namespace storm
-
-#endif /* STORM_UTILITY_IOUTILITY_H_ */
diff --git a/src/utility/Vector.h b/src/utility/Vector.h
index ab1fc461a..b3748d9ea 100644
--- a/src/utility/Vector.h
+++ b/src/utility/Vector.h
@@ -78,6 +78,18 @@ void subtractFromConstantOneVector(std::vector<T>* vector) {
 	}
 }
 
+template<class T>
+void addVectors(std::vector<T>& target, std::vector<T> const& summand) {
+	if (target.size() != target.size()) {
+		LOG4CPLUS_ERROR(logger, "Lengths of vectors does not match and makes operation impossible.");
+		throw storm::exceptions::InvalidArgumentException() << "Length of vectors does not match and makes operation impossible.";
+	}
+
+	for (uint_fast64_t i = 0; i < target.size(); ++i) {
+		target[i] += summand[i];
+	}
+}
+
 template<class T>
 void addVectors(std::vector<uint_fast64_t> const& states, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<T>& original, std::vector<T> const& summand) {
 	for (auto stateIt = states.cbegin(), stateIte = states.cend(); stateIt != stateIte; ++stateIt) {
diff --git a/test/functional/GmmxxDtmcPrctModelCheckerTest.cpp b/test/functional/GmmxxDtmcPrctModelCheckerTest.cpp
index 3050b3e49..059a7f295 100644
--- a/test/functional/GmmxxDtmcPrctModelCheckerTest.cpp
+++ b/test/functional/GmmxxDtmcPrctModelCheckerTest.cpp
@@ -17,7 +17,7 @@ TEST(GmmxxDtmcPrctModelCheckerTest, Die) {
 	ASSERT_EQ(dtmc->getNumberOfStates(), 13);
 	ASSERT_EQ(dtmc->getNumberOfTransitions(), 27);
 
-	storm::modelChecker::GmmxxDtmcPrctlModelChecker<double> mc(*dtmc);
+	storm::modelchecker::GmmxxDtmcPrctlModelChecker<double> mc(*dtmc);
 
 	storm::formula::Ap<double>* apFormula = new storm::formula::Ap<double>("one");
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(apFormula);
@@ -77,7 +77,7 @@ TEST(GmmxxDtmcPrctModelCheckerTest, Crowds) {
 	ASSERT_EQ(dtmc->getNumberOfStates(), 8607);
 	ASSERT_EQ(dtmc->getNumberOfTransitions(), 22460);
 
-	storm::modelChecker::GmmxxDtmcPrctlModelChecker<double> mc(*dtmc);
+	storm::modelchecker::GmmxxDtmcPrctlModelChecker<double> mc(*dtmc);
 
 	storm::formula::Ap<double>* apFormula = new storm::formula::Ap<double>("observe0Greater1");
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(apFormula);
@@ -125,7 +125,7 @@ TEST(GmmxxDtmcPrctModelCheckerTest, SynchronousLeader) {
 	ASSERT_EQ(dtmc->getNumberOfStates(), 12400);
 	ASSERT_EQ(dtmc->getNumberOfTransitions(), 28894);
 
-	storm::modelChecker::GmmxxDtmcPrctlModelChecker<double> mc(*dtmc);
+	storm::modelchecker::GmmxxDtmcPrctlModelChecker<double> mc(*dtmc);
 
 	storm::formula::Ap<double>* apFormula = new storm::formula::Ap<double>("elected");
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(apFormula);
diff --git a/test/functional/GmmxxMdpPrctModelCheckerTest.cpp b/test/functional/GmmxxMdpPrctModelCheckerTest.cpp
index 1fee30acd..ab6099d62 100644
--- a/test/functional/GmmxxMdpPrctModelCheckerTest.cpp
+++ b/test/functional/GmmxxMdpPrctModelCheckerTest.cpp
@@ -16,7 +16,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, Dice) {
 	ASSERT_EQ(mdp->getNumberOfStates(), 169);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), 436);
 
-	storm::modelChecker::GmmxxMdpPrctlModelChecker<double> mc(*mdp);
+	storm::modelchecker::GmmxxMdpPrctlModelChecker<double> mc(*mdp);
 
 	storm::formula::Ap<double>* apFormula = new storm::formula::Ap<double>("two");
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(apFormula);
@@ -112,7 +112,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, Dice) {
 
 	std::shared_ptr<storm::models::Mdp<double>> stateRewardMdp = stateRewardParser.getModel<storm::models::Mdp<double>>();
 
-	storm::modelChecker::GmmxxMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp);
+	storm::modelchecker::GmmxxMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp);
 
 	apFormula = new storm::formula::Ap<double>("done");
 	reachabilityRewardFormula = new storm::formula::ReachabilityReward<double>(apFormula);
@@ -142,7 +142,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, Dice) {
 
 	std::shared_ptr<storm::models::Mdp<double>> stateAndTransitionRewardMdp = stateAndTransitionRewardParser.getModel<storm::models::Mdp<double>>();
 
-	storm::modelChecker::GmmxxMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp);
+	storm::modelchecker::GmmxxMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp);
 
 	apFormula = new storm::formula::Ap<double>("done");
 	reachabilityRewardFormula = new storm::formula::ReachabilityReward<double>(apFormula);
@@ -178,7 +178,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, AsynchronousLeader) {
 	ASSERT_EQ(mdp->getNumberOfStates(), 3172);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), 7144);
 
-	storm::modelChecker::GmmxxMdpPrctlModelChecker<double> mc(*mdp);
+	storm::modelchecker::GmmxxMdpPrctlModelChecker<double> mc(*mdp);
 
 	storm::formula::Ap<double>* apFormula = new storm::formula::Ap<double>("elected");
 	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(apFormula);
diff --git a/test/parser/ParseDtmcTest.cpp b/test/parser/ParseDtmcTest.cpp
deleted file mode 100644
index d99c16f74..000000000
--- a/test/parser/ParseDtmcTest.cpp
+++ /dev/null
@@ -1,25 +0,0 @@
-/*
- * ParseDtmcTest.cpp
- *
- *  Created on: 03.12.2012
- *      Author: Thomas Heinemann
- */
-
-
-#include "gtest/gtest.h"
-#include "storm-config.h"
-#include "src/parser/DeterministicModelParser.h"
-#include "src/utility/IoUtility.h"
-
-TEST(ParseDtmcTest, parseAndOutput) {
-	storm::parser::DeterministicModelParser* dtmcParser = nullptr;
-	ASSERT_NO_THROW(dtmcParser = new storm::parser::DeterministicModelParser(
-			STORM_CPP_TESTS_BASE_PATH "/parser/tra_files/pctl_general_input_01.tra",
-			STORM_CPP_TESTS_BASE_PATH "/parser/lab_files/pctl_general_input_01.lab"));
-
-	ASSERT_NO_THROW(storm::utility::dtmcToDot(*(dtmcParser->getDtmc()), STORM_CPP_TESTS_BASE_PATH "/parser/output.dot"));
-
-	delete dtmcParser;
-}
-
-
diff --git a/test/parser/ParseMdpTest.cpp b/test/parser/ParseMdpTest.cpp
index 02302dc2a..b19335fdd 100644
--- a/test/parser/ParseMdpTest.cpp
+++ b/test/parser/ParseMdpTest.cpp
@@ -9,7 +9,6 @@
 #include "gtest/gtest.h"
 #include "storm-config.h"
 #include "src/parser/NondeterministicModelParser.h"
-#include "src/utility/IoUtility.h"
 
 TEST(ParseMdpTest, parseAndOutput) {
 	storm::parser::NondeterministicModelParser* mdpParser = nullptr;
diff --git a/test/parser/ReadLabFileTest.cpp b/test/parser/ReadLabFileTest.cpp
index 0e529104a..11ff3a71f 100644
--- a/test/parser/ReadLabFileTest.cpp
+++ b/test/parser/ReadLabFileTest.cpp
@@ -38,39 +38,42 @@ TEST(ReadLabFileTest, ParseTest) {
 		ASSERT_TRUE(labeling->containsAtomicProposition(smth));
 
 		//Testing whether all and only the correct nodes are labeled with "phi"
+		ASSERT_TRUE(labeling->stateHasAtomicProposition(phi,0));
 		ASSERT_TRUE(labeling->stateHasAtomicProposition(phi,1));
 		ASSERT_TRUE(labeling->stateHasAtomicProposition(phi,2));
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(phi,3));
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(phi,5));
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(phi,7));
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(phi,9));
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(phi,10));
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(phi,11));
 
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(phi,3));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(phi,4));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(phi,5));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(phi,6));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(phi,7));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(phi,8));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(phi,9));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(phi,10));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(phi,11));
 
 		//Testing whether all and only the correct nodes are labeled with "psi"
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(psi,6));
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(psi,7));
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(psi,8));
-
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,0));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,1));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,2));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,3));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,4));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,5));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,6));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,7));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,8));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,9));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,10));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(psi,11));
 
 		//Testing whether all and only the correct nodes are labeled with "smth"
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(smth,4));
-		ASSERT_TRUE(labeling->stateHasAtomicProposition(smth,5));
+		ASSERT_TRUE(labeling->stateHasAtomicProposition(smth,2));
 
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(smth,0));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(smth,1));
-		ASSERT_FALSE(labeling->stateHasAtomicProposition(smth,2));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(smth,3));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(smth,4));
+		ASSERT_FALSE(labeling->stateHasAtomicProposition(smth,5));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(smth,6));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(smth,7));
 		ASSERT_FALSE(labeling->stateHasAtomicProposition(smth,8));
diff --git a/test/parser/ReadTraFileTest.cpp b/test/parser/ReadTraFileTest.cpp
index cd75bc4a5..35ecfcbcc 100644
--- a/test/parser/ReadTraFileTest.cpp
+++ b/test/parser/ReadTraFileTest.cpp
@@ -12,8 +12,6 @@
 #include "src/exceptions/FileIoException.h"
 #include "src/exceptions/WrongFormatException.h"
 
-#include "src/utility/IoUtility.h"
-
 TEST(ReadTraFileTest, NonExistingFileTest) {
    //No matter what happens, please don't create a file with the name "nonExistingFile.not"! :-)
    ASSERT_THROW(storm::parser::DeterministicSparseTransitionParser(STORM_CPP_TESTS_BASE_PATH "/nonExistingFile.not"), storm::exceptions::FileIoException);
diff --git a/test/storage/SparseMatrixTest.cpp b/test/storage/SparseMatrixTest.cpp
index 96a1cea43..4ce2a7520 100644
--- a/test/storage/SparseMatrixTest.cpp
+++ b/test/storage/SparseMatrixTest.cpp
@@ -272,10 +272,6 @@ TEST(SparseMatrixTest, ConversionFromSparseEigen_RowMajor_SparseMatrixTest) {
 	ASSERT_NO_THROW(ssm = storm::adapters::StormAdapter::toStormSparseMatrix(esm));
 
 	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::ReadReady);
-	
-	const std::vector<uint_fast64_t> rowP = ssm->getRowIndications();
-	const std::vector<uint_fast64_t> colP = ssm->getColumnIndications();
-	const std::vector<int> valP = ssm->getStorage();
 
 	int target = -1;
 	for (auto &coeff: tripletList) {
diff --git a/test/storm-tests.cpp b/test/storm-tests.cpp
index 0fdcbb456..e5d78bb5e 100644
--- a/test/storm-tests.cpp
+++ b/test/storm-tests.cpp
@@ -37,7 +37,7 @@ void setUpLogging() {
 bool parseOptions(int const argc, char const * const argv[]) {
     storm::settings::Settings* s = nullptr;
     try {
-        storm::settings::Settings::registerModule<storm::modelChecker::GmmxxDtmcPrctlModelChecker<double>>();
+        storm::settings::Settings::registerModule<storm::modelchecker::GmmxxDtmcPrctlModelChecker<double>>();
         s = storm::settings::newInstance(argc, argv, nullptr, true);
     } catch (storm::exceptions::InvalidSettingsException& e) {
         std::cout << "Could not recover from settings error: " << e.what() << "." << std::endl;