Added native matrix-vector multiplication for our matrix format (as fast as gmm++). Fixed bug in bit vector. Fixed some issues in SCC decomposition. MDP model checkers now have the solving methods by default (native ones) and may override them with their own ones, if desired. Added some aux stuff, like vector helper methods.

src/modelchecker/MdpPrctlModelChecker.h

@@ -500,10 +500,85 @@ protected:
 	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 = 0;
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type>* summand = nullptr, uint_fast64_t repetitions = 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();
+
+		// Create vector for result of multiplication, which is reduced to the result vector after
+		// each multiplication.
+		std::vector<Type> multiplyResult(matrix.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);
+			}
+			if (this->minimumOperatorStack.top()) {
+				storm::utility::reduceVectorMin(multiplyResult, &vector, *nondeterministicChoiceIndices);
+			} else {
+				storm::utility::reduceVectorMax(multiplyResult, &vector, *nondeterministicChoiceIndices);
+			}
+		}
+	}
 
-	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const = 0;
+	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 {
+		// Get the settings object to customize solving.
+		storm::settings::Settings* s = storm::settings::instance();
+
+		// Get relevant user-defined settings for solving the equations.
+		double precision = s->get<double>("precision");
+		unsigned maxIterations = s->get<unsigned>("maxiter");
+		bool relative = s->get<bool>("relative");
+
+		// 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 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.
+		while (!converged && iterations < maxIterations) {
+			// Compute x' = A*x + b.
+			matrix.multiplyWithVector(*currentX, multiplyResult);
+
+			gmm::add(b, multiplyResult);
+
+			// Reduce the vector x' by applying min/max for all non-deterministic choices.
+			if (this->minimumOperatorStack.top()) {
+				storm::utility::reduceVectorMin(multiplyResult, newX, nondeterministicChoiceIndices);
+			} else {
+				storm::utility::reduceVectorMax(multiplyResult, newX, nondeterministicChoiceIndices);
+			}
 
+			// Determine whether the method converged.
+			converged = storm::utility::equalModuloPrecision(*currentX, *newX, precision, relative);
+
+			// Update environment variables.
+			swap = currentX;
+			currentX = newX;
+			newX = swap;
+			++iterations;
+		}
+
+		if (iterations % 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 " << iterations << " iterations.");
+		} else {
+			LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
+		}
+	}
 
 	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();

src/models/GraphTransitions.h

@@ -39,8 +39,8 @@ public:
 	 * @param forward If set to true, this objects will store the graph structure
 	 * of the backwards transition relation.
 	 */
-	GraphTransitions(std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix, bool forward)
-			: successorList(nullptr), stateIndications(nullptr), numberOfStates(transitionMatrix->getColumnCount()), numberOfTransitions(transitionMatrix->getNonZeroEntryCount()) {
+	GraphTransitions(storm::storage::SparseMatrix<T> const& transitionMatrix, bool forward)
+			: successorList(nullptr), stateIndications(nullptr), numberOfStates(transitionMatrix.getColumnCount()), numberOfTransitions(transitionMatrix.getNonZeroEntryCount()) {
 		if (forward) {
 			this->initializeForward(transitionMatrix);
 		} else {
@@ -48,12 +48,12 @@ public:
 		}
 	}
 
-	GraphTransitions(std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix, std::shared_ptr<std::vector<uint_fast64_t>> choiceIndices, bool forward)
-		: successorList(nullptr), stateIndications(nullptr), numberOfStates(transitionMatrix->getColumnCount()), numberOfTransitions(transitionMatrix->getNonZeroEntryCount()) {
+	GraphTransitions(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, bool forward)
+		: successorList(nullptr), stateIndications(nullptr), numberOfStates(transitionMatrix.getColumnCount()), numberOfTransitions(transitionMatrix.getNonZeroEntryCount()) {
 		if (forward) {
-			this->initializeForward(transitionMatrix, choiceIndices);
+			this->initializeForward(transitionMatrix, nondeterministicChoiceIndices);
 		} else {
-			this->initializeBackward(transitionMatrix, choiceIndices);
+			this->initializeBackward(transitionMatrix, nondeterministicChoiceIndices);
 		}
 	}
 
@@ -126,37 +126,37 @@ private:
 	 * Initializes this graph transitions object using the forward transition
 	 * relation given by means of a sparse matrix.
 	 */
-	void initializeForward(std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix) {
+	void initializeForward(storm::storage::SparseMatrix<T> const& transitionMatrix) {
 		this->successorList = new uint_fast64_t[numberOfTransitions];
 		this->stateIndications = new uint_fast64_t[numberOfStates + 1];
 
 		// 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);
+		std::copy(transitionMatrix.getRowIndications().begin(), transitionMatrix.getRowIndications().end(), this->stateIndications);
 
 		// 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.beginConstColumnIterator(i); rowIt != transitionMatrix.endConstColumnIterator(i); ++rowIt) {
 				this->successorList[currentNonZeroElement++] = *rowIt;
 			}
 		}
 	}
 
-	void initializeForward(std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix, std::shared_ptr<std::vector<uint_fast64_t>> choiceIndices) {
+	void initializeForward(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) {
 		this->successorList = new uint_fast64_t[numberOfTransitions];
 		this->stateIndications = new uint_fast64_t[numberOfStates + 1];
 
 		for (uint_fast64_t i = 0; i < numberOfStates; ++i) {
-			this->stateIndications[i] = transitionMatrix->getRowIndications().at(choiceIndices->at(i));
+			this->stateIndications[i] = transitionMatrix.getRowIndications().at(nondeterministicChoiceIndices[i]);
 		}
 		this->stateIndications[numberOfStates] = numberOfTransitions;
 
 		// 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 (uint_fast64_t j = choiceIndices->at(i); j < choiceIndices->at(i + 1); ++j) {
-				for (auto rowIt = transitionMatrix->beginConstColumnIterator(j); rowIt != transitionMatrix->endConstColumnIterator(j); ++rowIt) {
+			for (uint_fast64_t j = nondeterministicChoiceIndices[i]; j < nondeterministicChoiceIndices[i + 1]; ++j) {
+				for (auto rowIt = transitionMatrix.beginConstColumnIterator(j); rowIt != transitionMatrix.endConstColumnIterator(j); ++rowIt) {
 					this->successorList[currentNonZeroElement++] = *rowIt;
 				}
 			}
@@ -168,13 +168,13 @@ private:
 	 * relation, whose forward transition relation is given by means of a sparse
 	 * matrix.
 	 */
-	void initializeBackward(std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix) {
+	void initializeBackward(storm::storage::SparseMatrix<T> const& transitionMatrix) {
 		this->successorList = new uint_fast64_t[numberOfTransitions];
 		this->stateIndications = new uint_fast64_t[numberOfStates + 1]();
 
 		// 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.beginConstColumnIterator(i); rowIt != transitionMatrix.endConstColumnIterator(i); ++rowIt) {
 				this->stateIndications[*rowIt + 1]++;
 			}
 		}
@@ -197,7 +197,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.beginConstColumnIterator(i); rowIt != transitionMatrix.endConstColumnIterator(i); ++rowIt) {
 				this->successorList[nextIndicesList[*rowIt]++] = i;
 			}
 		}
@@ -206,14 +206,14 @@ private:
 		delete[] nextIndicesList;
 	}
 
-	void initializeBackward(std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix, std::shared_ptr<std::vector<uint_fast64_t>> choiceIndices) {
+	void initializeBackward(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) {
 		this->successorList = new uint_fast64_t[numberOfTransitions];
 		this->stateIndications = new uint_fast64_t[numberOfStates + 1]();
 
 		// 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 = choiceIndices->at(i); j < choiceIndices->at(i + 1); ++j) {
-				for (auto rowIt = transitionMatrix->beginConstColumnIterator(j); rowIt != transitionMatrix->endConstColumnIterator(j); ++rowIt) {
+			for (uint_fast64_t j = nondeterministicChoiceIndices[i]; j < nondeterministicChoiceIndices[i + 1]; ++j) {
+				for (auto rowIt = transitionMatrix.beginConstColumnIterator(j); rowIt != transitionMatrix.endConstColumnIterator(j); ++rowIt) {
 					this->stateIndications[*rowIt + 1]++;
 				}
 			}
@@ -237,8 +237,8 @@ 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 (uint_fast64_t j = (*choiceIndices)[i]; j < (*choiceIndices)[i + 1]; ++j) {
-				for (auto rowIt = transitionMatrix->beginConstColumnIterator(j); rowIt != transitionMatrix->endConstColumnIterator(j); ++rowIt) {
+			for (uint_fast64_t j = nondeterministicChoiceIndices[i]; j < nondeterministicChoiceIndices[i + 1]; ++j) {
+				for (auto rowIt = transitionMatrix.beginConstColumnIterator(j); rowIt != transitionMatrix.endConstColumnIterator(j); ++rowIt) {
 					this->successorList[nextIndicesList[*rowIt]++] = i;
 				}
 			}

src/storage/BitVector.h

@@ -547,7 +547,7 @@ private:
 		startingIndex >>= 6;
 		uint64_t* bucketPtr = this->bucketArray + startingIndex;
 
-		do {
+		while ((startingIndex << 6) < endIndex) {
 			// Compute the remaining bucket content by a right shift
 			// to the current bit.
 			uint_fast64_t remainingInBucket = *bucketPtr >> currentBitInByte;
@@ -570,7 +570,7 @@ private:
 
 			// Advance to the next bucket.
 			++startingIndex; ++bucketPtr; currentBitInByte = 0;
-		} while ((startingIndex << 6) < endIndex);
+		}
 		return endIndex;
 	}
 

src/storage/SparseMatrix.h

@@ -62,6 +62,44 @@ public:
 	 */
 	typedef const T* const constIterator;
 
+	class constRowsIterator {
+	public:
+		constRowsIterator(SparseMatrix<T> const& matrix) : matrix(matrix), offset(0) {
+			// Intentionally left empty.
+		}
+
+		constRowsIterator& operator++() {
+			++offset;
+			return *this;
+		}
+
+		bool operator==(constRowsIterator const& other) const {
+			return offset == other.offset;
+		}
+
+		bool operator!=(uint_fast64_t offset) const {
+			return this->offset != offset;
+		}
+
+		uint_fast64_t column() const {
+			return matrix.columnIndications[offset];
+		}
+
+		T const& value() const {
+			return matrix.valueStorage[offset];
+		}
+
+		void setOffset(uint_fast64_t offset) {
+			this->offset = offset;
+		}
+
+	private:
+		SparseMatrix<T> const& matrix;
+		uint_fast64_t offset;
+	};
+
+	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.
@@ -773,18 +811,44 @@ public:
 		return result;
 	}
 
-	T getRowVectorProduct(uint_fast64_t row, std::vector<T>& vector) {
-		T result = storm::utility::constGetZero<T>();;
-		auto valueIterator = valueStorage.begin() + rowIndications[row];
-		const auto valueIteratorEnd = valueStorage.begin() + rowIndications[row + 1];
-		auto columnIterator = columnIndications.begin() + rowIndications[row];
-		const auto columnIteratorEnd = columnIndications.begin() + rowIndications[row + 1];
-		for (; valueIterator != valueIteratorEnd; ++valueIterator, ++columnIterator) {
-			result += *valueIterator * vector[*columnIterator];
+	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 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;
+
+		std::cout << nondeterministicChoiceIndices << std::endl;
+		std::cout << rowCount << std::endl;
+
+		for (auto stateIt = states.cbegin(), stateIte = states.cend(); stateIt != stateIte; ++stateIt) {
+			std::cout << "stateIt " << *stateIt << std::endl;
+			std::cout << nondeterministicChoiceIndices[*stateIt] << std::endl;
+			std::cout << this->rowIndications[nondeterministicChoiceIndices[*stateIt]] << std::endl;
+			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);
+			}
+		}
+	}
+
 	/*!
 	 * Returns the size of the matrix in memory measured in bytes.
 	 * @return The size of the matrix in memory measured in bytes.
@@ -800,6 +864,10 @@ public:
 		return size;
 	}
 
+	constRowsIterator constRowsIteratorBegin() const {
+		return constRowsIterator(*this);
+	}
+
 	/*!
 	 * Returns an iterator to the columns of the non-zero entries of the given
 	 * row.
@@ -910,7 +978,7 @@ public:
 	 * to separate groups will be separated by a dashed line.
 	 * @return a (non-compressed) string representation of the matrix.
 	 */
-	std::string toString(std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices) const {
+	std::string toString(std::vector<uint_fast64_t> const* nondeterministicChoiceIndices) const {
 		std::stringstream result;
 		uint_fast64_t currentNondeterministicChoiceIndex = 0;
 

src/utility/GraphAnalyzer.h

@@ -71,7 +71,7 @@ public:
 		}
 
 		// Get the backwards transition relation from the model to ease the search.
-		storm::models::GraphTransitions<T> backwardTransitions(model.getTransitionMatrix(), false);
+		storm::models::GraphTransitions<T> backwardTransitions(*model.getTransitionMatrix(), false);
 
 		// Add all psi states as the already satisfy the condition.
 		*statesWithProbabilityGreater0 |= psiStates;
@@ -174,7 +174,7 @@ public:
 		}
 
 		// Get the backwards transition relation from the model to ease the search.
-		storm::models::GraphTransitions<T> backwardTransitions(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), false);
+		storm::models::GraphTransitions<T> backwardTransitions(*model.getTransitionMatrix(), *model.getNondeterministicChoiceIndices(), false);
 
 		// Add all psi states as the already satisfy the condition.
 		*statesWithProbability0 |= psiStates;
@@ -213,7 +213,7 @@ public:
 		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
 
 		// Get the backwards transition relation from the model to ease the search.
-		storm::models::GraphTransitions<T> backwardTransitions(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), false);
+		storm::models::GraphTransitions<T> backwardTransitions(*model.getTransitionMatrix(), *model.getNondeterministicChoiceIndices(), false);
 
 		storm::storage::BitVector* currentStates = new storm::storage::BitVector(model.getNumberOfStates(), true);
 
@@ -300,7 +300,7 @@ public:
 		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
 
 		// Get the backwards transition relation from the model to ease the search.
-		storm::models::GraphTransitions<T> backwardTransitions(transitionMatrix, nondeterministicChoiceIndices, false);
+		storm::models::GraphTransitions<T> backwardTransitions(*transitionMatrix, *nondeterministicChoiceIndices, false);
 
 		// Add all psi states as the already satisfy the condition.
 		*statesWithProbability0 |= psiStates;
@@ -360,7 +360,7 @@ public:
 		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
 
 		// Get the backwards transition relation from the model to ease the search.
-		storm::models::GraphTransitions<T> backwardTransitions(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), false);
+		storm::models::GraphTransitions<T> backwardTransitions(*model.getTransitionMatrix(), *model.getNondeterministicChoiceIndices(), false);
 
 		storm::storage::BitVector* currentStates = new storm::storage::BitVector(model.getNumberOfStates(), true);
 
@@ -417,14 +417,17 @@ public:
 	}
 
 	template <typename T>
-	static uint_fast64_t performSccDecomposition(storm::models::AbstractNondeterministicModel<T>& model, std::vector<std::vector<uint_fast64_t>*>* stronglyConnectedComponents) {
+	static uint_fast64_t performSccDecomposition(storm::storage::SparseMatrix<T> const& matrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<std::vector<uint_fast64_t>*>* stronglyConnectedComponents) {
 		LOG4CPLUS_INFO(logger, "Computing SCC decomposition.");
 
 		// Get the forward transition relation from the model to ease the search.
-		storm::models::GraphTransitions<T> forwardTransitions(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), true);
+		storm::models::GraphTransitions<T> forwardTransitions(matrix, nondeterministicChoiceIndices, true);
+
+		std::cout << matrix.toString(&nondeterministicChoiceIndices) << std::endl;
+		std::cout << forwardTransitions.toString() << std::endl;
 
 		// Perform the actual SCC decomposition based on the graph-transitions of the system.
-		uint_fast64_t result = performSccDecomposition(model.getNumberOfStates(), *model.getLabeledStates("init"), forwardTransitions, stronglyConnectedComponents);
+		uint_fast64_t result = performSccDecomposition(nondeterministicChoiceIndices.size(), forwardTransitions, stronglyConnectedComponents);
 
 		LOG4CPLUS_INFO(logger, "Done computing SCC decomposition.");
 		return result;
@@ -432,7 +435,7 @@ public:
 
 private:
 	template <typename T>
-	static uint_fast64_t performSccDecomposition(uint_fast64_t numberOfStates, storm::storage::BitVector const& initialStates, storm::models::GraphTransitions<T> const& forwardTransitions, std::vector<std::vector<uint_fast64_t>*>* stronglyConnectedComponents) {
+	static uint_fast64_t performSccDecomposition(uint_fast64_t numberOfStates, storm::models::GraphTransitions<T> const& forwardTransitions, std::vector<std::vector<uint_fast64_t>*>* stronglyConnectedComponents) {
 		std::vector<uint_fast64_t> tarjanStack;
 		tarjanStack.reserve(numberOfStates);
 		storm::storage::BitVector tarjanStackStates(numberOfStates);
@@ -442,7 +445,7 @@ private:
 		storm::storage::BitVector visitedStates(numberOfStates);
 
 		uint_fast64_t currentIndex = 0;
-		for (auto state : initialStates) {
+		for (uint_fast64_t state = 0; state < numberOfStates; ++state) {
 			if (!visitedStates.get(state)) {
 				performSccDecompositionHelper(state, currentIndex, stateIndices, lowlinks, tarjanStack, tarjanStackStates, visitedStates, forwardTransitions, stronglyConnectedComponents);
 			}

src/utility/Vector.h

@@ -78,6 +78,15 @@ void subtractFromConstantOneVector(std::vector<T>* vector) {
 	}
 }
 
+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) {
+		for (auto rowIt = nondeterministicChoiceIndices[*stateIt], rowIte = nondeterministicChoiceIndices[*stateIt + 1]; rowIt != rowIte; ++rowIt) {
+			original[rowIt] += summand[rowIt];
+		}
+	}
+}
+
 template<class T>
 void reduceVectorMin(std::vector<T> const& source, std::vector<T>* target, std::vector<uint_fast64_t> const& filter) {
 	uint_fast64_t currentSourceRow = 0;