diff --git a/examples/mdp/scc/scc.pctl b/examples/mdp/scc/scc.pctl
new file mode 100644
index 000000000..393670a26
--- /dev/null
+++ b/examples/mdp/scc/scc.pctl
@@ -0,0 +1,2 @@
+Pmin=? [ F a ]
+Pmax=? [ F a ]
\ No newline at end of file
diff --git a/src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h b/src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h
index 750d068c8..7093f4cb3 100644
--- a/src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h
@@ -9,6 +9,7 @@
 #define STORM_MODELCHECKER_PRCTL_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_H_
 
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
+#include "src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.h"
 #include "src/exceptions/InvalidPropertyException.h"
 #include <cmath>
 
@@ -29,7 +30,7 @@ public:
 	 *
 	 * @param model The MDP to be checked.
 	 */
-	explicit TopologicalValueIterationMdpPrctlModelChecker(storm::models::Mdp<Type> const& model) : SparseMdpPrctlModelChecker<Type>(model) {
+	explicit TopologicalValueIterationMdpPrctlModelChecker(storm::models::Mdp<Type> const& model) : SparseMdpPrctlModelChecker<Type>(model, std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<Type>>(new storm::solver::TopologicalValueIterationNondeterministicLinearEquationSolver<Type>())) {
 		// Intentionally left empty.
 	}
 
@@ -46,100 +47,6 @@ public:
 	 * Virtual destructor. Needs to be virtual, because this class has virtual methods.
 	 */
 	virtual ~TopologicalValueIterationMdpPrctlModelChecker() { }
-
-private:
-	/*!
-	 * Solves the given equation system under the given parameters using the power method.
-	 *
-	 * @param A The matrix A specifying the coefficients of the 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.
-	 */
-	void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
-		// Get the settings object to customize solving.
-		storm::settings::Settings* s = storm::settings::Settings::getInstance();
-
-		// Get relevant user-defined settings for solving the equations.
-		double precision = s->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
-		uint_fast64_t maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger();
-		bool relative = s->getOptionByLongName("relative").getArgument(0).getValueAsBoolean();
-
-		// Now, we need to determine the SCCs of the MDP and a topological sort.
-        std::vector<std::vector<uint_fast64_t>> stronglyConnectedComponents = storm::utility::graph::performSccDecomposition(matrix, nondeterministicChoiceIndices, stronglyConnectedComponents, stronglyConnectedComponentsDependencyGraph);
-        storm::storage::SparseMatrix<T> stronglyConnectedComponentsDependencyGraph = this->getModel().extractSccDependencyGraph(stronglyConnectedComponents);
-		std::vector<uint_fast64_t> topologicalSort = storm::utility::graph::getTopologicalSort(stronglyConnectedComponentsDependencyGraph);
-
-		// Set up the environment for the power method.
-		std::vector<Type> multiplyResult(matrix.getRowCount());
-		std::vector<Type>* currentX = &x;
-		std::vector<Type>* newX = new std::vector<Type>(x.size());
-		std::vector<Type>* swap = nullptr;
-		uint_fast64_t currentMaxLocalIterations = 0;
-		uint_fast64_t localIterations = 0;
-		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) {
-				// Compute x' = A*x + b.
-				matrix.multiplyWithVector(scc, nondeterministicChoiceIndices, *currentX, multiplyResult);
-				storm::utility::addVectors(scc, nondeterministicChoiceIndices, multiplyResult, b);
-
-				// Reduce the vector x' by applying min/max for all non-deterministic choices.
-				if (this->minimumOperatorStack.top()) {
-					storm::utility::reduceVectorMin(multiplyResult, newX, scc, nondeterministicChoiceIndices);
-				} else {
-					storm::utility::reduceVectorMax(multiplyResult, newX, scc, nondeterministicChoiceIndices);
-				}
-
-				// 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);
-
-				// Update environment variables.
-				swap = currentX;
-				currentX = newX;
-				newX = swap;
-				++localIterations;
-				++globalIterations;
-			}
-
-			// 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;
-		} else {
-			delete newX;
-		}
-
-		// Check if the solver converged and issue a warning otherwise.
-		if (converged) {
-			LOG4CPLUS_INFO(logger, "Iterative solver converged after " << currentMaxLocalIterations << " iterations.");
-		} else {
-			LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
-		}
-	}
 };
 
 } // namespace prctl
diff --git a/src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.cpp b/src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.cpp
index e075e2ee6..6ff0310ba 100644
--- a/src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.cpp
+++ b/src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.cpp
@@ -66,16 +66,24 @@ namespace storm {
 
 			// 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.
+			int counter = 0;
+			std::cout << "Solving Equation System using the TopologicalValueIterationNon..." << std::endl;
 			for (auto sccIndexIt = topologicalSort.begin(); sccIndexIt != topologicalSort.end() && converged; ++sccIndexIt) {
-				std::vector<uint_fast64_t> const& scc = sccDecomposition[*sccIndexIt];
+				storm::storage::StateBlock const& scc = sccDecomposition[*sccIndexIt];
+
+				std::cout << "SCC " << counter << " contains:" << std::endl;
+				for (auto sccIt = scc.cbegin(); sccIt != scc.cend(); ++sccIt) {
+					std::cout << *sccIt << ", ";
+				}
+				std::cout << std::endl;
 
 				// For the current SCC, we need to perform value iteration until convergence.
 				localIterations = 0;
 				converged = false;
 				while (!converged && localIterations < maximalNumberOfIterations) {
 					// Compute x' = A*x + b.
-					A.multiplyWithVector(scc, nondeterministicChoiceIndices, *currentX, multiplyResult);
-					storm::utility::addVectors(scc, nondeterministicChoiceIndices, multiplyResult, b);
+					//A.multiplyWithVector(scc, nondeterministicChoiceIndices, *currentX, multiplyResult);
+					//storm::utility::addVectors(scc, nondeterministicChoiceIndices, multiplyResult, b);
 
 					/*
 					Versus:
@@ -85,17 +93,17 @@ namespace storm {
 
 					// Reduce the vector x' by applying min/max for all non-deterministic choices.
 					if (minimize) {
-						storm::utility::reduceVectorMin(*multiplyResult, *newX, scc, nondeterministicChoiceIndices);
+						//storm::utility::reduceVectorMin(*multiplyResult, *newX, scc, nondeterministicChoiceIndices);
 					}
 					else {
-						storm::utility::reduceVectorMax(*multiplyResult, *newX, scc, nondeterministicChoiceIndices);
+						//storm::utility::reduceVectorMax(*multiplyResult, *newX, scc, nondeterministicChoiceIndices);
 					}
 
 					// 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);
+					converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, precision, relative);
 
 					// Update environment variables.
 					swap = currentX;
@@ -120,7 +128,7 @@ namespace storm {
 			}
 			
 			if (!xMemoryProvided) {
-				delete copyX;
+				delete newX;
 			}
 
 			if (!multiplyResultMemoryProvided) {
diff --git a/src/utility/graph.h b/src/utility/graph.h
index 8cb7c044e..9459014bd 100644
--- a/src/utility/graph.h
+++ b/src/utility/graph.h
@@ -569,7 +569,7 @@ namespace storm {
                     LOG4CPLUS_ERROR(logger, "Provided matrix is required to be square.");
                     throw storm::exceptions::InvalidArgumentException() << "Provided matrix is required to be square.";
                 }
-                
+
                 uint_fast64_t numberOfStates = matrix.getRowCount();
                 
                 // Prepare the result. This relies on the matrix being square.
@@ -598,12 +598,12 @@ namespace storm {
                             
                         recursionStepBackward:
                             for (; successorIterator != matrix.end(currentState); ++successorIterator) {
-                                if (!visitedStates.get(successorIterator.first)) {
+                                if (!visitedStates.get(successorIterator->first)) {
                                     // Put unvisited successor on top of our recursion stack and remember that.
-                                    recursionStack.push_back(successorIterator.first);
+									recursionStack.push_back(successorIterator->first);
                                     
                                     // Also, put initial value for iterator on corresponding recursion stack.
-                                    iteratorRecursionStack.push_back(matrix.begin(successorIterator.first));
+									iteratorRecursionStack.push_back(matrix.begin(successorIterator->first));
                                     
                                     goto recursionStepForward;
                                 }
diff --git a/test/functional/modelchecker/TopologicalValueIterationMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/TopologicalValueIterationMdpPrctlModelCheckerTest.cpp
index 7ccb05369..114de1b47 100644
--- a/test/functional/modelchecker/TopologicalValueIterationMdpPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/TopologicalValueIterationMdpPrctlModelCheckerTest.cpp
@@ -9,8 +9,9 @@
 
 TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, Dice) {
 	storm::settings::Settings* s = storm::settings::Settings::getInstance();
-	storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", "", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.trans.rew");
-    
+	//storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", "", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.trans.rew");
+	storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/mdp/scc/scc.tra", STORM_CPP_BASE_PATH "/examples/mdp/scc/scc.lab", "");
+
 	ASSERT_EQ(parser.getType(), storm::models::MDP);
     
 	std::shared_ptr<storm::models::Mdp<double>> mdp = parser.getModel<storm::models::Mdp<double>>();
@@ -152,7 +153,7 @@ TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, Dice) {
 }
 
 TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, AsynchronousLeader) {
-	storm::settings::Settings* s = storm::settings::Settings::getInstance();
+	/*storm::settings::Settings* s = storm::settings::Settings::getInstance();
 	storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/mdp/asynchronous_leader/leader4.tra", STORM_CPP_BASE_PATH "/examples/mdp/asynchronous_leader/leader4.lab", "", STORM_CPP_BASE_PATH "/examples/mdp/asynchronous_leader/leader4.trans.rew");
 
 	ASSERT_EQ(parser.getType(), storm::models::MDP);
@@ -220,5 +221,5 @@ TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, AsynchronousLeader) {
 	result = mc.checkNoBoundOperator(*rewardFormula);;
 
 	ASSERT_LT(std::abs(result[0] - 4.285689611), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-	delete rewardFormula;
+	delete rewardFormula;*/
 }