diff --git a/src/modelchecker/MdpPrctlModelChecker.h b/src/modelchecker/MdpPrctlModelChecker.h
index 8b48f11a3..acc75b2fb 100644
--- a/src/modelchecker/MdpPrctlModelChecker.h
+++ b/src/modelchecker/MdpPrctlModelChecker.h
@@ -545,6 +545,7 @@ private:
 		while (!converged && iterations < maxIterations) {
 			// Compute x' = A*x + b.
 			matrix.multiplyWithVector(*currentX, multiplyResult);
+			// matrix.multiplyAddAndReduceInPlace(nondeterministicChoiceIndices, *currentX, b, this->minimumOperatorStack.top());
 
 			gmm::add(b, multiplyResult);
 
@@ -558,11 +559,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 (iterations % 2 == 1) {
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index afdb0ea5b..5259a8052 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -830,6 +830,17 @@ public:
 		}
 	}
 
+	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;
@@ -843,6 +854,37 @@ public:
 		}
 	}
 
+	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;
+				}
+			}
+		}
+	}
+
 	/*!
 	 * Returns the size of the matrix in memory measured in bytes.
 	 * @return The size of the matrix in memory measured in bytes.