LRA finally working for ctmcs

Former-commit-id: 699e4714a4

src/modelchecker/csl/SparseCtmcCslModelChecker.cpp

@@ -12,6 +12,8 @@
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 
+#include "src/storage/StronglyConnectedComponentDecomposition.h"
+
 #include "src/exceptions/InvalidStateException.h"
 #include "src/exceptions/InvalidPropertyException.h"
 #include "src/exceptions/NotImplementedException.h"
@@ -476,25 +478,193 @@ namespace storm {
             std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
 
-            // Since we use the uniformized matrix to compute the steady state probabilities, we need to build it first.
-            ValueType uniformizationRate = 0;
-            for (auto const& rate : this->getModel().getExitRateVector()) {
-                uniformizationRate = std::max(uniformizationRate, rate);
+            storm::storage::SparseMatrix<ValueType> probabilityMatrix = this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector());
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(computeLongRunAverageHelper(this->getModel(), probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector(), qualitative, *linearEquationSolverFactory)));
+        }
+        
+        template<typename ValueType>
+        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(storm::models::sparse::DeterministicModel<ValueType> const& model, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            // If there are no goal states, we avoid the computation and directly return zero.
+            uint_fast64_t numOfStates = transitionMatrix.getRowCount();
+            if (psiStates.empty()) {
+                return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>());
             }
-            uniformizationRate *= 1.02;
-            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
             
-//            storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), uniformizationRate, this->getModel().getExitRateVector());
-
-//            storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeGeneratorMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector());
-            storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector());
+            // Likewise, if all bits are set, we can avoid the computation.
+            if (psiStates.full()) {
+                return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>());
+            }
+            
+            // Start by decomposing the DTMC into its BSCCs.
+            storm::storage::StronglyConnectedComponentDecomposition<double> bsccDecomposition(model, false, true);
+            
+            // Get some data members for convenience.
+            ValueType one = storm::utility::one<ValueType>();
+            ValueType zero = storm::utility::zero<ValueType>();
             
-            std::vector<ValueType> result = SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(this->getModel(), uniformizedMatrix, subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory);
+            // First we check which states are in BSCCs.
+            storm::storage::BitVector statesInBsccs(numOfStates);
+            storm::storage::BitVector firstStatesInBsccs(numOfStates);
             
+            std::vector<uint_fast64_t> stateToBsccIndexMap(transitionMatrix.getColumnCount());
             
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>());
+            for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
+                storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
+                
+                // Gather information for later use.
+                bool first = true;
+                for (auto const& state : bscc) {
+                    statesInBsccs.set(state);
+                    if (first) {
+                        firstStatesInBsccs.set(state);
+                    }
+                    stateToBsccIndexMap[state] = currentBsccIndex;
+                    first = false;
+                }
+            }
+            
+            firstStatesInBsccs = firstStatesInBsccs % statesInBsccs;
+            storm::storage::BitVector statesNotInBsccs = ~statesInBsccs;
+            
+            // Then we construct an equation system that yields the steady state probabilities for all states in BSCCs.
+            storm::storage::SparseMatrix<ValueType> bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true);
+            
+            // Since in the fix point equation, we need to multiply the vector from the left, we convert this to a
+            // multiplication from the right by transposing the system.
+            bsccEquationSystem = bsccEquationSystem.transpose(false, true);
+            
+            // Create an auxiliary structure that makes it easy to look up the indices within the set of BSCC states.
+            uint_fast64_t lastIndex = 0;
+            uint_fast64_t currentNumberOfSetBits = 0;
+            std::vector<uint_fast64_t> indexInStatesInBsccs;
+            indexInStatesInBsccs.reserve(transitionMatrix.getRowCount());
+            for (auto index : statesInBsccs) {
+                while (lastIndex <= index) {
+                    indexInStatesInBsccs.push_back(currentNumberOfSetBits);
+                    ++lastIndex;
+                }
+                ++currentNumberOfSetBits;
+            }
+            
+            // Now build the final equation system matrix.
+            storm::storage::SparseMatrixBuilder<ValueType> builder;
+            uint_fast64_t currentBsccIndex = 0;
+            for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) {
+                // If the current row is the first one belonging to a BSCC, we substitute it by the constraint that the
+                // values for states of this BSCC must sum to one.
+                if (firstStatesInBsccs.get(row)) {
+                    storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
+                    
+                    for (auto const& state : bscc) {
+                        builder.addNextValue(row, indexInStatesInBsccs[state], one);
+                    }
+                    ++currentBsccIndex;
+                } else {
+                    // Otherwise, we copy the row, and subtract 1 from the diagonal.
+                    for (auto& entry : bsccEquationSystem.getRow(row)) {
+                        if (entry.getColumn() == row) {
+                            builder.addNextValue(row, entry.getColumn(), entry.getValue() - one);
+                        } else {
+                            builder.addNextValue(row, entry.getColumn(), entry.getValue());
+                        }
+                    }
+                }
+            }
+            
+            bsccEquationSystem = builder.build();
+            
+            std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero);
+            storm::utility::vector::setVectorValues(bsccEquationSystemRightSide, firstStatesInBsccs, one);
+            
+            // As an initial guess, we take the uniform distribution over all states of the containing BSCC.
+            std::vector<ValueType> bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), zero);
+            for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
+                storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
+                
+                for (auto const& state : bscc) {
+                    bsccEquationSystemSolution[indexInStatesInBsccs[state]] = one / bscc.size();
+                }
+            }
+            
+            {
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(bsccEquationSystem);
+                solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide);
+            }
+            
+            // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour.
+            if (exitRateVector != nullptr) {
+                std::vector<ValueType> bsccTotalValue(bsccDecomposition.size(), zero);
+                std::size_t i = 0;
+                for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) {
+                    bsccTotalValue[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i] * (one / (*exitRateVector)[*stateIter]);
+                }
+                
+                i = 0;
+                for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) {
+                    bsccEquationSystemSolution[i] = (bsccEquationSystemSolution[i] * (one / (*exitRateVector)[*stateIter])) / bsccTotalValue[stateToBsccIndexMap[*stateIter]];
+                }
+            }
+            
+            // Calculate LRA Value for each BSCC from steady state distribution in BSCCs.
+            std::vector<ValueType> bsccLra(bsccDecomposition.size(), zero);
+            size_t i = 0;
+            for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) {
+                if (psiStates.get(*stateIter)) {
+                    bsccLra[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i];
+                }
+            }
+            
+            std::vector<ValueType> rewardSolution;
+            if (!statesNotInBsccs.empty()) {
+                // Calculate LRA for states not in bsccs as expected reachability rewards.
+                // Target states are states in bsccs, transition reward is the lra of the bscc for each transition into a
+                // bscc and 0 otherwise. This corresponds to the sum of LRAs in BSCC weighted by the reachability probability
+                // of the BSCC.
+                
+                std::vector<ValueType> rewardRightSide;
+                rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits());
+                
+                for (auto state : statesNotInBsccs) {
+                    ValueType reward = zero;
+                    for (auto entry : transitionMatrix.getRow(state)) {
+                        if (statesInBsccs.get(entry.getColumn())) {
+                            reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]];
+                        }
+                    }
+                    rewardRightSide.push_back(reward);
+                }
+                
+                storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true);
+                rewardEquationSystemMatrix.convertToEquationSystem();
+                
+                rewardSolution = std::vector<ValueType>(rewardEquationSystemMatrix.getColumnCount(), one);
+                
+                {
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(rewardEquationSystemMatrix);
+                    solver->solveEquationSystem(rewardSolution, rewardRightSide);
+                }
+            }
+            
+            // Fill the result vector.
+            std::vector<ValueType> result(numOfStates);
+            auto rewardSolutionIter = rewardSolution.begin();
+            for (size_t state = 0; state < numOfStates; ++state) {
+                if (statesInBsccs.get(state)) {
+                    // Assign the value of the bscc the state is in.
+                    result[state] = bsccLra[stateToBsccIndexMap[state]];
+                } else {
+                    STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution.");
+                    // Take the value from the reward computation. Since the n-th state not in any bscc is the n-th
+                    // entry in rewardSolution we can just take the next value from the iterator.
+                    result[state] = *rewardSolutionIter;
+                    ++rewardSolutionIter;
+                }
+            }
+            
+            return result;
         }
         
+        
         // Explicitly instantiate the model checker.
         template class SparseCtmcCslModelChecker<double>;
         

src/modelchecker/csl/SparseCtmcCslModelChecker.h

@@ -13,10 +13,14 @@ namespace storm {
         template<storm::dd::DdType DdType, typename ValueType>
         class HybridCtmcCslModelChecker;
         
+        template<typename ValueType>
+        class SparseDtmcPrctlModelChecker;
+        
         template<class ValueType>
         class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker<ValueType> {
         public:
             friend class HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, ValueType>;
+            friend class SparseDtmcPrctlModelChecker<ValueType>;
             
             explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model);
             explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
@@ -40,7 +44,7 @@ namespace storm {
             static std::vector<ValueType> computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             std::vector<ValueType> computeInstantaneousRewardsHelper(double timeBound) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(double timeBound) const;
-
+            static std::vector<ValueType> computeLongRunAverageHelper(storm::models::sparse::DeterministicModel<ValueType> const& model, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             /*!
              * Computes the matrix representing the transitions of the uniformized CTMC.
              *

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp

@@ -8,11 +8,11 @@
 #include "src/utility/graph.h"
 #include "src/utility/solver.h"
 
+#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 
 #include "src/exceptions/InvalidStateException.h"
-#include "src/storage/StronglyConnectedComponentDecomposition.h"
 
 #include "src/exceptions/InvalidPropertyException.h"
 
@@ -302,280 +302,13 @@ namespace storm {
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative)));
         }
-        
-		template<typename ValueType>
-        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(storm::models::sparse::DeterministicModel<ValueType> const& model, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-			// If there are no goal states, we avoid the computation and directly return zero.
-			uint_fast64_t numOfStates = transitionMatrix.getRowCount();
-			if (psiStates.empty()) {
-				return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>());
-			}
-
-			// Likewise, if all bits are set, we can avoid the computation.
-			if (psiStates.full()) {
-				return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>());
-			}
-
-			// Start by decomposing the DTMC into its BSCCs.
-			storm::storage::StronglyConnectedComponentDecomposition<double> bsccDecomposition(model, false, true);
-
-			// Get some data members for convenience.
-			ValueType one = storm::utility::one<ValueType>();
-			ValueType zero = storm::utility::zero<ValueType>();
-
-			// First we check which states are in BSCCs. We use this later to speed up reachability analysis
-			storm::storage::BitVector statesInBsccs(numOfStates);
-            storm::storage::BitVector statesInBsccsWithoutFirst(numOfStates);
-    
-			std::vector<uint_fast64_t> stateToBsccIndexMap(transitionMatrix.getColumnCount());
-
-			for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
-				storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
-
-				// Gather information for later use.
-                bool first = true;
-				for (auto const& state : bscc) {
-					statesInBsccs.set(state);
-                    if (!first) {
-                        statesInBsccsWithoutFirst.set(state);
-                    }
-					stateToBsccIndexMap[state] = currentBsccIndex;
-                    first = false;
-				}
-			}
-
-			storm::storage::BitVector statesNotInBsccs = ~statesInBsccs;
-
-            std::cout << transitionMatrix << std::endl;
-            
-			// Calculate steady state distribution for all BSCCs by calculating an eigenvector for the eigenvalue 1 of
-            // the transposed transition matrix for the BSCCs.
-			storm::storage::SparseMatrix<ValueType> bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true);
-            
-            std::cout << bsccEquationSystem << std::endl;
-            
-            bsccEquationSystem = bsccEquationSystem.transpose(false, true);
-
-            std::cout << bsccEquationSystem << std::endl;
-            
-			// Subtract identity matrix.
-			for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) {
-				for (auto& entry : bsccEquationSystem.getRow(row)) {
-					if (entry.getColumn() == row) {
-						entry.setValue(entry.getValue() - one);
-					}
-				}
-			}
-            
-            std::cout << bsccEquationSystem << std::endl;
-            
-            std::cout << statesInBsccsWithoutFirst << " // " << statesInBsccs << std::endl;
-            bsccEquationSystem = bsccEquationSystem.getSubmatrix(false, statesInBsccsWithoutFirst, statesInBsccs, false);
-            
-            std::cout << bsccEquationSystem << std::endl;
-
-            // For each BSCC, add a row to the matrix that expresses that the sum over all entries in this BSCC needs to be one.
-            storm::storage::SparseMatrixBuilder<ValueType> builder(std::move(bsccEquationSystem));
-            typename storm::storage::SparseMatrixBuilder<ValueType>::index_type row = builder.getLastRow() + 1;
-            
-            for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
-                storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
-                
-                for (auto const& state : bscc) {
-                    builder.addNextValue(row, state, one);
-                }
-                ++row;
-            }
-            bsccEquationSystem = builder.build();
-            
-            std::cout << bsccEquationSystem << std::endl;
-            
-			std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero);
-            bsccEquationSystemRightSide.back() = one;
-			std::vector<ValueType> bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one);
-			{
-                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(bsccEquationSystem);
-				solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide);
-			}
-
-            ValueType sum = zero;
-            for (auto const& elem : bsccEquationSystemSolution) {
-                std::cout << "sol " << elem << std::endl;
-                sum += elem;
-            }
-            std::cout << "sum: " << sum << std::endl;
-            std::cout << "in " << bsccDecomposition.size() << "bsccs" << std::endl;
-            
-			// Calculate LRA Value for each BSCC from steady state distribution in BSCCs.
-			// We have to scale the results, as the probabilities for each BSCC have to sum up to one, which they don't
-            // necessarily do in the solution of the equation system.
-			std::vector<ValueType> bsccLra(bsccDecomposition.size(), zero);
-			std::vector<ValueType> bsccTotalValue(bsccDecomposition.size(), zero);
-			size_t i = 0;
-			for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) {
-				if (psiStates.get(*stateIter)) {
-					bsccLra[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i];
-				}
-				bsccTotalValue[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i];
-			}
-			for (i = 0; i < bsccLra.size(); ++i) {
-				bsccLra[i] /= bsccTotalValue[i];
-			}
-            
-            std::vector<ValueType> rewardSolution;
-            if (!statesNotInBsccs.empty()) {
-                // Calculate LRA for states not in bsccs as expected reachability rewards.
-                // Target states are states in bsccs, transition reward is the lra of the bscc for each transition into a
-                // bscc and 0 otherwise. This corresponds to the sum of LRAs in BSCC weighted by the reachability probability
-                // of the BSCC.
-                
-                std::vector<ValueType> rewardRightSide;
-                rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits());
-                
-                for (auto state : statesNotInBsccs) {
-                    ValueType reward = zero;
-                    for (auto entry : transitionMatrix.getRow(state)) {
-                        if (statesInBsccs.get(entry.getColumn())) {
-                            reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]];
-                        }
-                    }
-                    rewardRightSide.push_back(reward);
-                }
-                
-                storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true);
-                rewardEquationSystemMatrix.convertToEquationSystem();
-                
-                rewardSolution = std::vector<ValueType>(rewardEquationSystemMatrix.getColumnCount(), one);
-                
-                {
-                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(rewardEquationSystemMatrix);
-                    solver->solveEquationSystem(rewardSolution, rewardRightSide);
-                }
-            }
-            
-            // Fill the result vector.
-            std::vector<ValueType> result(numOfStates);
-            auto rewardSolutionIter = rewardSolution.begin();
-            for (size_t state = 0; state < numOfStates; ++state) {
-                if (statesInBsccs.get(state)) {
-                    // Assign the value of the bscc the state is in.
-                    result[state] = bsccLra[stateToBsccIndexMap[state]];
-                } else {
-                    STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution.");
-                    // Take the value from the reward computation. Since the n-th state not in any bscc is the n-th
-                    // entry in rewardSolution we can just take the next value from the iterator.
-                    result[state] = *rewardSolutionIter;
-                    ++rewardSolutionIter;
-                }
-            }
-
-			return result;
-
-			//old implementeation
-
-			//now we calculate the long run average for each BSCC in isolation and compute the weighted contribution of the BSCC to the LRA value of all states
-			//for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
-			//	storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
-
-			//	storm::storage::BitVector statesInThisBscc(numOfStates);
-			//	for (auto const& state : bscc) {
-			//		statesInThisBscc.set(state);
-			//	}
-
-			//	//ValueType lraForThisBscc = this->computeLraForBSCC(transitionMatrix, psiStates, bscc);
-			//	ValueType lraForThisBscc = bsccLra[currentBsccIndex];
-
-			//	//the LRA value of a BSCC contributes to the LRA value of a state with the probability of reaching that BSCC from that state
-			//	//thus we add Prob(Eventually statesInThisBscc) * lraForThisBscc to the result vector
-
-			//	//the reachability probabilities will be zero in other BSCCs, thus we can set the left operand of the until formula to statesNotInBsccs as an optimization
-			//	std::vector<ValueType> reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBscc, false);
-			//	
-			//	storm::utility::vector::scaleVectorInPlace<ValueType>(reachProbThisBscc, lraForThisBscc);
-			//	storm::utility::vector::addVectorsInPlace<ValueType>(result, reachProbThisBscc);
-			//}
-
-			//return result;
-		}
 
 		template<typename ValueType>
 		std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
 			std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
 			ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
 
-			return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeLongRunAverageHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory)));
-		}
-
-
-		template<typename ValueType>
-		ValueType SparseDtmcPrctlModelChecker<ValueType>::computeLraForBSCC(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
-			//if size is one we can immediately derive the result
-			if (bscc.size() == 1){
-				if (psiStates.get(*(bscc.begin()))) {
-					return storm::utility::one<ValueType>();
-				} else{
-					return storm::utility::zero<ValueType>();
-				}
-			}
-
-			storm::storage::BitVector subsystem = storm::storage::BitVector(transitionMatrix.getRowCount());
-			subsystem.set(bscc.begin(), bscc.end());
-
-			//we now have to solve ((P')^t - I) * x = 0, where P' is the submatrix of transitionMatrix,
-			// ^t means transose, and I is the identity matrix.
-			
-			storm::storage::SparseMatrix<ValueType> subsystemMatrix = transitionMatrix.getSubmatrix(false, subsystem, subsystem, true);
-			subsystemMatrix = subsystemMatrix.transpose();
-
-			// subtract 1 on the diagonal and at the same time add a row with all ones to enforce that the result of the equation system is unique
-			storm::storage::SparseMatrixBuilder<ValueType> equationSystemBuilder(subsystemMatrix.getRowCount() + 1, subsystemMatrix.getColumnCount(), subsystemMatrix.getEntryCount() + subsystemMatrix.getColumnCount());
-			ValueType one = storm::utility::one<ValueType>();
-			ValueType zero = storm::utility::zero<ValueType>();
-			bool foundDiagonalElement = false;
-			for (uint_fast64_t row = 0; row < subsystemMatrix.getRowCount(); ++row) {
-				for (auto& entry : subsystemMatrix.getRowGroup(row)) {
-					if (entry.getColumn() == row) {
-						equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue() - one);
-						foundDiagonalElement = true;
-					} else {
-						equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue());
-					}
-				}
-
-				// Throw an exception if a row did not have an element on the diagonal.
-				STORM_LOG_THROW(foundDiagonalElement, storm::exceptions::InvalidOperationException, "Internal Error, no diagonal entry found.");
-			}
-			for (uint_fast64_t column = 0; column + 1 < subsystemMatrix.getColumnCount(); ++column) {
-				equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), column, one);
-			}
-			equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), subsystemMatrix.getColumnCount() - 1, zero);
-			subsystemMatrix = equationSystemBuilder.build();
-
-			// create x and b for the equation system. setting the last entry of b to one enforces the sum over the unique solution vector is one
-			std::vector<ValueType> steadyStateDist(subsystemMatrix.getRowCount(), zero);
-			std::vector<ValueType> b(subsystemMatrix.getRowCount(), zero);
-			b[subsystemMatrix.getRowCount() - 1] = one;
-
-
-			{
-				auto solver = linearEquationSolverFactory.create(subsystemMatrix);
-				solver->solveEquationSystem(steadyStateDist, b);
-			}
-
-			//remove the last entry of the vector as it was just there to enforce the unique solution
-			steadyStateDist.pop_back();
-			
-			//calculate the fraction we spend in psi states on the long run
-			std::vector<ValueType> steadyStateForPsiStates(transitionMatrix.getRowCount() - 1, zero);
-			storm::utility::vector::setVectorValues(steadyStateForPsiStates, psiStates & subsystem, steadyStateDist);
-
-			ValueType result = zero;
-
-			for (auto value : steadyStateForPsiStates) {
-				result += value;
-			}
-
-			return result;
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseCtmcCslModelChecker<ValueType>::computeLongRunAverageHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), nullptr, qualitative, *linearEquationSolverFactory)));
 		}
 
         template<typename ValueType>

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h

@@ -43,8 +43,6 @@ namespace storm {
             static std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
 			static std::vector<ValueType> computeLongRunAverageHelper(storm::models::sparse::DeterministicModel<ValueType> const& model, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
 
-			static ValueType computeLraForBSCC(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-
             // An object that is used for retrieving linear equation solvers.
             std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };

src/settings/modules/GmmxxEquationSolverSettings.cpp

@@ -7,7 +7,7 @@ namespace storm {
         namespace modules {
             
             const std::string GmmxxEquationSolverSettings::moduleName = "gmm++";
-            const std::string GmmxxEquationSolverSettings::techniqueOptionName = "tech";
+            const std::string GmmxxEquationSolverSettings::techniqueOptionName = "method";
             const std::string GmmxxEquationSolverSettings::preconditionOptionName = "precond";
             const std::string GmmxxEquationSolverSettings::restartOptionName = "restart";
             const std::string GmmxxEquationSolverSettings::maximalIterationsOptionName = "maxiter";
@@ -32,38 +32,38 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build());
             }
             
-            bool GmmxxEquationSolverSettings::isLinearEquationSystemTechniqueSet() const {
+            bool GmmxxEquationSolverSettings::isLinearEquationSystemMethodSet() const {
                 return this->getOption(techniqueOptionName).getHasOptionBeenSet();
             }
             
-            GmmxxEquationSolverSettings::LinearEquationTechnique GmmxxEquationSolverSettings::getLinearEquationSystemTechnique() const {
+            GmmxxEquationSolverSettings::LinearEquationMethod GmmxxEquationSolverSettings::getLinearEquationSystemMethod() const {
                 std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString();
                 if (linearEquationSystemTechniqueAsString == "bicgstab") {
-                    return GmmxxEquationSolverSettings::LinearEquationTechnique::Bicgstab;
+                    return GmmxxEquationSolverSettings::LinearEquationMethod::Bicgstab;
                 } else if (linearEquationSystemTechniqueAsString == "qmr") {
-                    return GmmxxEquationSolverSettings::LinearEquationTechnique::Qmr;
+                    return GmmxxEquationSolverSettings::LinearEquationMethod::Qmr;
                 } else if (linearEquationSystemTechniqueAsString == "gmres") {
-                    return GmmxxEquationSolverSettings::LinearEquationTechnique::Gmres;
+                    return GmmxxEquationSolverSettings::LinearEquationMethod::Gmres;
                 } else if (linearEquationSystemTechniqueAsString == "jacobi") {
-                    return GmmxxEquationSolverSettings::LinearEquationTechnique::Jacobi;
+                    return GmmxxEquationSolverSettings::LinearEquationMethod::Jacobi;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected.");
             }
             
-            bool GmmxxEquationSolverSettings::isPreconditioningTechniqueSet() const {
+            bool GmmxxEquationSolverSettings::isPreconditioningMethodSet() const {
                 return this->getOption(preconditionOptionName).getHasOptionBeenSet();
             }
             
-            GmmxxEquationSolverSettings::PreconditioningTechnique GmmxxEquationSolverSettings::getPreconditioningTechnique() const {
-                std::string preconditioningTechniqueAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString();
-                if (preconditioningTechniqueAsString == "ilu") {
-                    return GmmxxEquationSolverSettings::PreconditioningTechnique::Ilu;
-                } else if (preconditioningTechniqueAsString == "diagonal") {
-                    return GmmxxEquationSolverSettings::PreconditioningTechnique::Diagonal;
-                } else if (preconditioningTechniqueAsString == "none") {
-                    return GmmxxEquationSolverSettings::PreconditioningTechnique::None;
+            GmmxxEquationSolverSettings::PreconditioningMethod GmmxxEquationSolverSettings::getPreconditioningMethod() const {
+                std::string PreconditioningMethodAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString();
+                if (PreconditioningMethodAsString == "ilu") {
+                    return GmmxxEquationSolverSettings::PreconditioningMethod::Ilu;
+                } else if (PreconditioningMethodAsString == "diagonal") {
+                    return GmmxxEquationSolverSettings::PreconditioningMethod::Diagonal;
+                } else if (PreconditioningMethodAsString == "none") {
+                    return GmmxxEquationSolverSettings::PreconditioningMethod::None;
                 }
-                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << preconditioningTechniqueAsString << "' selected.");
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << PreconditioningMethodAsString << "' selected.");
             }
             
             bool GmmxxEquationSolverSettings::isRestartIterationCountSet() const {
@@ -99,7 +99,7 @@ namespace storm {
             }
             
             bool GmmxxEquationSolverSettings::check() const {
-                bool optionsSet = isLinearEquationSystemTechniqueSet() || isPreconditioningTechniqueSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet();
+                bool optionsSet = isLinearEquationSystemMethodSet() || isPreconditioningMethodSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet();
                 
                 STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx || !optionsSet, "gmm++ is not selected as the equation solver, so setting options for gmm++ has no effect.");
                 

src/settings/modules/GmmxxEquationSolverSettings.h

@@ -12,11 +12,11 @@ namespace storm {
              */
             class GmmxxEquationSolverSettings : public ModuleSettings {
             public:
-                // An enumeration of all available techniques for solving linear equations.
-                enum class LinearEquationTechnique { Bicgstab, Qmr, Gmres, Jacobi };
+                // An enumeration of all available methods for solving linear equations.
+                enum class LinearEquationMethod { Bicgstab, Qmr, Gmres, Jacobi };
 
-                // An enumeration of all available preconditioning techniques.
-                enum class PreconditioningTechnique { Ilu, Diagonal, None };
+                // An enumeration of all available preconditioning methods.
+                enum class PreconditioningMethod { Ilu, Diagonal, None };
                 
                 // An enumeration of all available convergence criteria.
                 enum class ConvergenceCriterion { Absolute, Relative };
@@ -29,32 +29,32 @@ namespace storm {
                 GmmxxEquationSolverSettings(storm::settings::SettingsManager& settingsManager);
                 
                 /*!
-                 * Retrieves whether the linear equation system technique has been set.
+                 * Retrieves whether the linear equation system method has been set.
                  *
-                 * @return True iff the linear equation system technique has been set.
+                 * @return True iff the linear equation system method has been set.
                  */
-                bool isLinearEquationSystemTechniqueSet() const;
+                bool isLinearEquationSystemMethodSet() const;
                 
                 /*!
-                 * Retrieves the technique that is to be used for solving systems of linear equations.
+                 * Retrieves the method that is to be used for solving systems of linear equations.
                  *
-                 * @return The technique to use.
+                 * @return The method to use.
                  */
-                LinearEquationTechnique getLinearEquationSystemTechnique() const;
+                LinearEquationMethod getLinearEquationSystemMethod() const;
                 
                 /*!
-                 * Retrieves whether the preconditioning technique has been set.
+                 * Retrieves whether the preconditioning method has been set.
                  *
-                 * @return True iff the preconditioning technique has been set.
+                 * @return True iff the preconditioning method has been set.
                  */
-                bool isPreconditioningTechniqueSet() const;
+                bool isPreconditioningMethodSet() const;
                 
                 /*!
-                 * Retrieves the technique that is to be used for preconditioning solving systems of linear equations.
+                 * Retrieves the method that is to be used for preconditioning solving systems of linear equations.
                  *
-                 * @return The technique to use.
+                 * @return The method to use.
                  */
-                PreconditioningTechnique getPreconditioningTechnique() const;
+                PreconditioningMethod getPreconditioningMethod() const;
                 
                 /*!
                  * Retrieves whether the restart iteration count has been set.
@@ -64,7 +64,7 @@ namespace storm {
                 bool isRestartIterationCountSet() const;
                 
                 /*!
-                 * Retrieves the number of iterations after which restarted techniques are to be restarted.
+                 * Retrieves the number of iterations after which restarted methods are to be restarted.
                  *
                  * @return The number of iterations after which to restart.
                  */

src/settings/modules/NativeEquationSolverSettings.cpp

@@ -7,20 +7,23 @@ namespace storm {
         namespace modules {
             
             const std::string NativeEquationSolverSettings::moduleName = "native";
-            const std::string NativeEquationSolverSettings::techniqueOptionName = "tech";
+            const std::string NativeEquationSolverSettings::techniqueOptionName = "method";
+            const std::string NativeEquationSolverSettings::omegaOptionName = "soromega";
             const std::string NativeEquationSolverSettings::maximalIterationsOptionName = "maxiter";
             const std::string NativeEquationSolverSettings::maximalIterationsOptionShortName = "i";
             const std::string NativeEquationSolverSettings::precisionOptionName = "precision";
             const std::string NativeEquationSolverSettings::absoluteOptionName = "absolute";
             
             NativeEquationSolverSettings::NativeEquationSolverSettings(storm::settings::SettingsManager& settingsManager) : ModuleSettings(settingsManager, moduleName) {
-                std::vector<std::string> methods = { "jacobi" };
+                std::vector<std::string> methods = { "jacobi", "gaussseidel", "sor" };
                 this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the native engine. Available are: { jacobi }.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("jacobi").build()).build());
                 
                 this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build());
                 
                 this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
                 
+                this->addOption(storm::settings::OptionBuilder(moduleName, omegaOptionName, false, "The omega used for SOR.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The value of the SOR parameter.").setDefaultValueDouble(0.9).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+                
                 this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build());
             }
             
@@ -28,10 +31,14 @@ namespace storm {
                 return this->getOption(techniqueOptionName).getHasOptionBeenSet();
             }
             
-            NativeEquationSolverSettings::LinearEquationTechnique NativeEquationSolverSettings::getLinearEquationSystemTechnique() const {
+            NativeEquationSolverSettings::LinearEquationMethod NativeEquationSolverSettings::getLinearEquationSystemMethod() const {
                 std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString();
                 if (linearEquationSystemTechniqueAsString == "jacobi") {
-                    return NativeEquationSolverSettings::LinearEquationTechnique::Jacobi;
+                    return NativeEquationSolverSettings::LinearEquationMethod::Jacobi;
+                } else if (linearEquationSystemTechniqueAsString == "gaussseidel") {
+                    return NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel;
+                } else if (linearEquationSystemTechniqueAsString == "sor") {
+                    return NativeEquationSolverSettings::LinearEquationMethod::SOR;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected.");
             }
@@ -51,6 +58,10 @@ namespace storm {
             double NativeEquationSolverSettings::getPrecision() const {
                 return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble();
             }
+
+            double NativeEquationSolverSettings::getOmega() const {
+                return this->getOption(omegaOptionName).getArgumentByName("value").getValueAsDouble();
+            }
             
             bool NativeEquationSolverSettings::isConvergenceCriterionSet() const {
                 return this->getOption(absoluteOptionName).getHasOptionBeenSet();

src/settings/modules/NativeEquationSolverSettings.h

@@ -12,8 +12,8 @@ namespace storm {
              */
             class NativeEquationSolverSettings : public ModuleSettings {
             public:
-                // An enumeration of all available techniques for solving linear equations.
-                enum class LinearEquationTechnique { Jacobi };
+                // An enumeration of all available methods for solving linear equations.
+                enum class LinearEquationMethod { Jacobi, GaussSeidel, SOR };
                 
                 // An enumeration of all available convergence criteria.
                 enum class ConvergenceCriterion { Absolute, Relative };
@@ -33,11 +33,11 @@ namespace storm {
                 bool isLinearEquationSystemTechniqueSet() const;
                 
                 /*!
-                 * Retrieves the technique that is to be used for solving systems of linear equations.
+                 * Retrieves the method that is to be used for solving systems of linear equations.
                  *
-                 * @return The technique to use.
+                 * @return The method to use.
                  */
-                LinearEquationTechnique getLinearEquationSystemTechnique() const;
+                LinearEquationMethod getLinearEquationSystemMethod() const;
                 
                 /*!
                  * Retrieves whether the maximal iteration count has been set.
@@ -67,6 +67,13 @@ namespace storm {
                  */
                 double getPrecision() const;
                 
+                /*!
+                 * Retrieves the value of omega to be used for the SOR method.
+                 *
+                 * @return The value of omega to be used.
+                 */
+                double getOmega() const;
+                
                 /*!
                  * Retrieves whether the convergence criterion has been set.
                  *
@@ -89,6 +96,7 @@ namespace storm {
             private:
                 // Define the string names of the options as constants.
                 static const std::string techniqueOptionName;
+                static const std::string omegaOptionName;
                 static const std::string maximalIterationsOptionName;
                 static const std::string maximalIterationsOptionShortName;
                 static const std::string precisionOptionName;

src/solver/GmmxxLinearEquationSolver.cpp

@@ -31,24 +31,24 @@ namespace storm {
             restart = settings.getRestartIterationCount();
             
             // Determine the method to be used.
-            storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique methodAsSetting = settings.getLinearEquationSystemTechnique();
-            if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Bicgstab) {
+            storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod methodAsSetting = settings.getLinearEquationSystemMethod();
+            if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Bicgstab) {
                 method = SolutionMethod::Bicgstab;
-            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Qmr) {
+            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Qmr) {
                 method = SolutionMethod::Qmr;
-            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Gmres) {
+            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Gmres) {
                 method = SolutionMethod::Gmres;
-            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Jacobi) {
+            } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Jacobi) {
                 method = SolutionMethod::Jacobi;
             }
             
             // Check which preconditioner to use.
-            storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique preconditionAsSetting = settings.getPreconditioningTechnique();
-            if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::Ilu) {
+            storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod preconditionAsSetting = settings.getPreconditioningMethod();
+            if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::Ilu) {
                 preconditioner = Preconditioner::Ilu;
-            } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::Diagonal) {
+            } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::Diagonal) {
                 preconditioner = Preconditioner::Diagonal;
-            } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::None) {
+            } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::None) {
                 preconditioner = Preconditioner::None;
             }
         }
@@ -149,12 +149,6 @@ namespace storm {
         uint_fast64_t GmmxxLinearEquationSolver<ValueType>::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             // Get a Jacobi decomposition of the matrix A.
             std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
-            std::cout << "LU" << std::endl;
-            std::cout << jacobiDecomposition.first << std::endl;
-            for (auto const& elem : jacobiDecomposition.second) {
-                std::cout << elem << std::endl;
-            }
-            std::cout << "----" << std::endl;
             
             // Convert the LU matrix to gmm++'s format.
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.first));

src/solver/NativeLinearEquationSolver.cpp

@@ -15,7 +15,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : A(A) {
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method) : A(A), method(method) {
             // Get the settings object to customize linear solving.
             storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
             
@@ -23,61 +23,94 @@ namespace storm {
             maximalNumberOfIterations = settings.getMaximalIterationCount();
             precision = settings.getPrecision();
             relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
-            
-            // Determine the method to be used.
-            storm::settings::modules::NativeEquationSolverSettings::LinearEquationTechnique methodAsSetting = settings.getLinearEquationSystemTechnique();
-            if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationTechnique::Jacobi) {
-                method = SolutionMethod::Jacobi;
-            }
         }
                 
         template<typename ValueType>
         void NativeLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
-            // Get a Jacobi decomposition of the matrix A.
-            std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
-            
-            // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
-            bool multiplyResultProvided = true;
-            std::vector<ValueType>* nextX = multiplyResult;
-            if (nextX == nullptr) {
-                nextX = new std::vector<ValueType>(x.size());
-                multiplyResultProvided = false;
-            }
-            std::vector<ValueType> const* copyX = nextX;
-            std::vector<ValueType>* currentX = &x;
-            
-            // Target vector for precision calculation.
-            std::vector<ValueType> tmpX(x.size());
-            
-            // Set up additional environment variables.
-            uint_fast64_t iterationCount = 0;
-            bool converged = false;
-            
-            while (!converged && iterationCount < maximalNumberOfIterations) {
-                // Compute D^-1 * (b - LU * x) and store result in nextX.
-                jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
-                storm::utility::vector::subtractVectors(b, tmpX, tmpX);
-                storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
+            if (method == SolutionMethod::SOR || method == SolutionMethod::GaussSeidel) {
+                // Define the omega used for SOR.
+                ValueType omega = method == SolutionMethod::SOR ? storm::settings::nativeEquationSolverSettings().getOmega() : storm::utility::one<ValueType>();
                 
-                // Swap the two pointers as a preparation for the next iteration.
-                std::swap(nextX, currentX);
+                // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
+                bool tmpXProvided = true;
+                std::vector<ValueType>* tmpX = multiplyResult;
+                if (multiplyResult == nullptr) {
+                    tmpX = new std::vector<ValueType>(x);
+                    tmpXProvided = false;
+                } else {
+                    *tmpX = x;
+                }
                 
-                // Now check if the process already converged within our precision.
-                converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(precision), relative);
+                // Set up additional environment variables.
+                uint_fast64_t iterationCount = 0;
+                bool converged = false;
                 
-                // Increase iteration count so we can abort if convergence is too slow.
-                ++iterationCount;
-            }
-            
-            // 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 input parameter x.
-            if (currentX == copyX) {
-                std::swap(x, *currentX);
-            }
-            
-            // If the vector for the temporary multiplication result was not provided, we need to delete it.
-            if (!multiplyResultProvided) {
-                delete copyX;
+                while (!converged && iterationCount < maximalNumberOfIterations) {
+                    A.performSuccessiveOverRelaxationStep(omega, x, b);
+                    
+                    // Now check if the process already converged within our precision.
+                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(x, *tmpX, static_cast<ValueType>(precision), relative);
+                    
+                    // If we did not yet converge, we need to copy the contents of x to *tmpX.
+                    if (!converged) {
+                        *tmpX = x;
+                    }
+                    
+                    // Increase iteration count so we can abort if convergence is too slow.
+                    ++iterationCount;
+                }
+                
+                // If the vector for the temporary multiplication result was not provided, we need to delete it.
+                if (!tmpXProvided) {
+                    delete tmpX;
+                }
+            } else {
+                // Get a Jacobi decomposition of the matrix A.
+                std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+                
+                // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
+                bool multiplyResultProvided = true;
+                std::vector<ValueType>* nextX = multiplyResult;
+                if (nextX == nullptr) {
+                    nextX = new std::vector<ValueType>(x.size());
+                    multiplyResultProvided = false;
+                }
+                std::vector<ValueType> const* copyX = nextX;
+                std::vector<ValueType>* currentX = &x;
+                
+                // Target vector for precision calculation.
+                std::vector<ValueType> tmpX(x.size());
+                
+                // Set up additional environment variables.
+                uint_fast64_t iterationCount = 0;
+                bool converged = false;
+                
+                while (!converged && iterationCount < maximalNumberOfIterations) {
+                    // Compute D^-1 * (b - LU * x) and store result in nextX.
+                    jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
+                    storm::utility::vector::subtractVectors(b, tmpX, tmpX);
+                    storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
+                    
+                    // Swap the two pointers as a preparation for the next iteration.
+                    std::swap(nextX, currentX);
+                    
+                    // Now check if the process already converged within our precision.
+                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(precision), relative);
+                    
+                    // Increase iteration count so we can abort if convergence is too slow.
+                    ++iterationCount;
+                }
+                
+                // 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 input parameter x.
+                if (currentX == copyX) {
+                    std::swap(x, *currentX);
+                }
+                
+                // If the vector for the temporary multiplication result was not provided, we need to delete it.
+                if (!multiplyResultProvided) {
+                    delete copyX;
+                }
             }
         }
         

src/solver/NativeLinearEquationSolver.h

@@ -14,15 +14,16 @@ namespace storm {
         public:
             // An enumeration specifying the available solution methods.
             enum class SolutionMethod {
-                Jacobi
+                Jacobi, GaussSeidel, SOR
             };
             
             /*!
              * Constructs a linear equation solver with parameters being set according to the settings object.
              *
              * @param A The matrix defining the coefficients of the linear equation system.
+             * @param method The method to use for solving linear equations.
              */
-            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
+            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method = SolutionMethod::Jacobi);
 
             /*!
              * Constructs a linear equation solver with the given parameters.

src/storage/SparseMatrix.cpp

@@ -50,7 +50,7 @@ namespace storm {
         void MatrixEntry<IndexType, ValueType>::setValue(ValueType const& value) {
             this->entry.second = value;
         }
-
+        
         template<typename IndexType, typename ValueType>
         std::pair<IndexType, ValueType> const& MatrixEntry<IndexType, ValueType>::getColumnValuePair() const {
             return this->entry;
@@ -102,7 +102,7 @@ namespace storm {
             if (row < lastRow) {
                 throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrixBuilder::addNextValue: adding an element in row " << row << ", but an element in row " << lastRow << " has already been added.";
             }
-
+            
             // Check that we did not move backwards wrt. to column.
             if (row == lastRow && column < lastColumn) {
                 throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrixBuilder::addNextValue: adding an element in column " << column << " in row " << row << ", but an element in column " << lastColumn << " has already been added in that row.";
@@ -117,14 +117,14 @@ namespace storm {
                 
                 lastRow = row;
             }
-                        
+            
             lastColumn = column;
             
             // Finally, set the element and increase the current size.
             columnsAndValues.emplace_back(column, value);
             highestColumn = std::max(highestColumn, column);
             ++currentEntryCount;
-
+            
             // In case we did not expect this value, we throw an exception.
             if (forceInitialDimensions) {
                 STORM_LOG_THROW(!initialRowCountSet || lastRow < initialRowCount, storm::exceptions::OutOfRangeException, "Cannot insert value at illegal row " << lastRow << ".");
@@ -159,14 +159,14 @@ namespace storm {
             // as now the indices of row i are always between rowIndications[i] and rowIndications[i + 1], also for
             // the first and last row.
             rowIndications.push_back(currentEntryCount);
-
+            
             uint_fast64_t columnCount = highestColumn + 1;
             if (initialColumnCountSet && forceInitialDimensions) {
                 STORM_LOG_THROW(columnCount <= initialColumnCount, storm::exceptions::InvalidStateException, "Expected not more than " << initialColumnCount << " columns, but got " << columnCount << ".");
                 columnCount = std::max(columnCount, initialColumnCount);
             }
             columnCount = std::max(columnCount, overriddenColumnCount);
-
+            
             uint_fast64_t entryCount = currentEntryCount;
             if (initialEntryCountSet && forceInitialDimensions) {
                 STORM_LOG_THROW(entryCount == initialEntryCount, storm::exceptions::InvalidStateException, "Expected " << initialEntryCount << " entries, but got " << entryCount << ".");
@@ -184,12 +184,12 @@ namespace storm {
                     rowGroupCount = std::max(rowGroupCount, initialRowGroupCount);
                 }
                 rowGroupCount = std::max(rowGroupCount, overriddenRowGroupCount);
-
+                
                 for (index_type i = currentRowGroup; i <= rowGroupCount; ++i) {
                     rowGroupIndices.push_back(rowCount);
                 }
             }
-
+            
             return SparseMatrix<ValueType>(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), hasCustomRowGrouping);
         }
         
@@ -270,14 +270,14 @@ namespace storm {
         
         template<typename ValueType>
         SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, ValueType>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(rowGroupIndices) {
-			this->updateNonzeroEntryCount();
+            this->updateNonzeroEntryCount();
         }
         
         template<typename ValueType>
         SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, ValueType>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(std::move(rowGroupIndices)) {
-			this->updateNonzeroEntryCount();
+            this->updateNonzeroEntryCount();
         }
-                
+        
         template<typename ValueType>
         SparseMatrix<ValueType>& SparseMatrix<ValueType>::operator=(SparseMatrix<ValueType> const& other) {
             // Only perform assignment if source and target are not the same.
@@ -321,7 +321,7 @@ namespace storm {
             }
             
             bool equalityResult = true;
-
+            
             equalityResult &= this->getRowCount() == other.getRowCount();
             if (!equalityResult) {
                 return false;
@@ -379,35 +379,35 @@ namespace storm {
             return entryCount;
         }
         
-		template<typename T>
-		uint_fast64_t SparseMatrix<T>::getRowGroupEntryCount(uint_fast64_t const group) const {
-			uint_fast64_t result = 0;
-			for (uint_fast64_t row = this->getRowGroupIndices()[group]; row < this->getRowGroupIndices()[group + 1]; ++row) {
-				result += (this->rowIndications[row + 1] - this->rowIndications[row]);
-			}
-			return result;
-		}
-
+        template<typename T>
+        uint_fast64_t SparseMatrix<T>::getRowGroupEntryCount(uint_fast64_t const group) const {
+            uint_fast64_t result = 0;
+            for (uint_fast64_t row = this->getRowGroupIndices()[group]; row < this->getRowGroupIndices()[group + 1]; ++row) {
+                result += (this->rowIndications[row + 1] - this->rowIndications[row]);
+            }
+            return result;
+        }
+        
         template<typename ValueType>
         typename SparseMatrix<ValueType>::index_type SparseMatrix<ValueType>::getNonzeroEntryCount() const {
             return nonzeroEntryCount;
-		}
-
-		template<typename ValueType>
-		void SparseMatrix<ValueType>::updateNonzeroEntryCount() const {
-			//SparseMatrix<ValueType>* self = const_cast<SparseMatrix<ValueType>*>(this);
-			this->nonzeroEntryCount = 0;
-			for (auto const& element : *this) {
-				if (element.getValue() != storm::utility::zero<ValueType>()) {
-					++this->nonzeroEntryCount;
-				}
-			}
-		}
-
-		template<typename ValueType>
-		void SparseMatrix<ValueType>::updateNonzeroEntryCount(std::make_signed<index_type>::type difference) {
-			this->nonzeroEntryCount += difference;
-		}
+        }
+        
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::updateNonzeroEntryCount() const {
+            //SparseMatrix<ValueType>* self = const_cast<SparseMatrix<ValueType>*>(this);
+            this->nonzeroEntryCount = 0;
+            for (auto const& element : *this) {
+                if (element.getValue() != storm::utility::zero<ValueType>()) {
+                    ++this->nonzeroEntryCount;
+                }
+            }
+        }
+        
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::updateNonzeroEntryCount(std::make_signed<index_type>::type difference) {
+            this->nonzeroEntryCount += difference;
+        }
         
         template<typename ValueType>
         typename SparseMatrix<ValueType>::index_type SparseMatrix<ValueType>::getRowGroupCount() const {
@@ -423,7 +423,7 @@ namespace storm {
         std::vector<typename SparseMatrix<ValueType>::index_type> const& SparseMatrix<ValueType>::getRowGroupIndices() const {
             return rowGroupIndices;
         }
-
+        
         template<typename ValueType>
         void SparseMatrix<ValueType>::makeRowsAbsorbing(storm::storage::BitVector const& rows) {
             for (auto row : rows) {
@@ -457,7 +457,7 @@ namespace storm {
             columnValuePtr->setValue(storm::utility::one<ValueType>());
             ++columnValuePtr;
             for (; columnValuePtr != columnValuePtrEnd; ++columnValuePtr) {
-				++this->nonzeroEntryCount;
+                ++this->nonzeroEntryCount;
                 columnValuePtr->setColumn(0);
                 columnValuePtr->setValue(storm::utility::zero<ValueType>());
             }
@@ -510,33 +510,42 @@ namespace storm {
                 return getSubmatrix(rowConstraint, columnConstraint, fakeRowGroupIndices, insertDiagonalElements);
             }
         }
-                
+        
         template<typename ValueType>
         SparseMatrix<ValueType> SparseMatrix<ValueType>::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector<index_type> const& rowGroupIndices, bool insertDiagonalEntries) const {
             uint_fast64_t submatrixColumnCount = columnConstraint.getNumberOfSetBits();
             
             // Start by creating a temporary vector that stores for each index whose bit is set to true the number of
             // bits that were set before that particular index.
-            std::vector<index_type> bitsSetBeforeIndex;
-            bitsSetBeforeIndex.reserve(columnCount);
+            std::vector<index_type> columnBitsSetBeforeIndex;
+            columnBitsSetBeforeIndex.reserve(columnCount);
             
             // Compute the information to fill this vector.
             index_type lastIndex = 0;
             index_type currentNumberOfSetBits = 0;
-            
-            // If we are requested to add missing diagonal entries, we need to make sure the corresponding rows are also
-            // taken.
-//            storm::storage::BitVector columnBitCountConstraint = columnConstraint;
-//            if (insertDiagonalEntries) {
-//                columnBitCountConstraint |= rowGroupConstraint;
-//            }
             for (auto index : columnConstraint) {
                 while (lastIndex <= index) {
-                    bitsSetBeforeIndex.push_back(currentNumberOfSetBits);
+                    columnBitsSetBeforeIndex.push_back(currentNumberOfSetBits);
                     ++lastIndex;
                 }
                 ++currentNumberOfSetBits;
             }
+            std::vector<index_type>* rowBitsSetBeforeIndex;
+            if (&rowGroupConstraint == &columnConstraint) {
+                rowBitsSetBeforeIndex = &columnBitsSetBeforeIndex;
+            } else {
+                rowBitsSetBeforeIndex = new std::vector<index_type>(rowCount);
+                
+                lastIndex = 0;
+                currentNumberOfSetBits = 0;
+                for (auto index : rowGroupConstraint) {
+                    while (lastIndex <= index) {
+                        rowBitsSetBeforeIndex->push_back(currentNumberOfSetBits);
+                        ++lastIndex;
+                    }
+                    ++currentNumberOfSetBits;
+                }
+            }
             
             // Then, we need to determine the number of entries and the number of rows of the submatrix.
             index_type subEntries = 0;
@@ -551,7 +560,7 @@ namespace storm {
                         if (columnConstraint.get(it->getColumn())) {
                             ++subEntries;
                             
-                            if (index == it->getColumn()) {
+                            if (columnBitsSetBeforeIndex[it->getColumn()] == (*rowBitsSetBeforeIndex)[index]) {
                                 foundDiagonalElement = true;
                             }
                         }
@@ -580,13 +589,13 @@ namespace storm {
                     
                     for (const_iterator it = this->begin(i), ite = this->end(i); it != ite; ++it) {
                         if (columnConstraint.get(it->getColumn())) {
-                            if (index == it->getColumn()) {
+                            if (columnBitsSetBeforeIndex[it->getColumn()] == (*rowBitsSetBeforeIndex)[index]) {
                                 insertedDiagonalElement = true;
-                            } else if (insertDiagonalEntries && !insertedDiagonalElement && it->getColumn() > index) {
-                                matrixBuilder.addNextValue(rowCount, rowCount, storm::utility::zero<ValueType>());
+                            } else if (insertDiagonalEntries && !insertedDiagonalElement && columnBitsSetBeforeIndex[it->getColumn()] > (*rowBitsSetBeforeIndex)[index]) {
+                                matrixBuilder.addNextValue(rowCount, rowGroupCount, storm::utility::zero<ValueType>());
                                 insertedDiagonalElement = true;
                             }
-                            matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[it->getColumn()], it->getValue());
+                            matrixBuilder.addNextValue(rowCount, columnBitsSetBeforeIndex[it->getColumn()], it->getValue());
                         }
                     }
                     if (insertDiagonalEntries && !insertedDiagonalElement && rowGroupCount < submatrixColumnCount) {
@@ -597,6 +606,11 @@ namespace storm {
                 ++rowGroupCount;
             }
             
+            // If the constraints were not the same, we have allocated some additional memory we need to free now.
+            if (&rowGroupConstraint != &columnConstraint) {
+                delete rowBitsSetBeforeIndex;
+            }
+            
             return matrixBuilder.build();
         }
         
@@ -655,12 +669,12 @@ namespace storm {
         SparseMatrix<ValueType> SparseMatrix<ValueType>::transpose(bool joinGroups, bool keepZeros) const {
             index_type rowCount = this->getColumnCount();
             index_type columnCount = joinGroups ? this->getRowGroupCount() : this->getRowCount();
-			if (keepZeros) {
-				index_type entryCount = this->getEntryCount();
-			} else {
-				this->updateNonzeroEntryCount();
-				index_type entryCount = this->getNonzeroEntryCount();
-			}
+            if (keepZeros) {
+                index_type entryCount = this->getEntryCount();
+            } else {
+                this->updateNonzeroEntryCount();
+                index_type entryCount = this->getNonzeroEntryCount();
+            }
             
             std::vector<index_type> rowIndications(rowCount + 1);
             std::vector<MatrixEntry<index_type, ValueType>> columnsAndValues(entryCount);
@@ -668,7 +682,7 @@ namespace storm {
             // First, we need to count how many entries each column has.
             for (index_type group = 0; group < columnCount; ++group) {
                 for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) {
-					if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
+                    if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
                         ++rowIndications[transition.getColumn() + 1];
                     }
                 }
@@ -687,7 +701,7 @@ namespace storm {
             // Now we are ready to actually fill in the values of the transposed matrix.
             for (index_type group = 0; group < columnCount; ++group) {
                 for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) {
-					if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
+                    if (transition.getValue() != storm::utility::zero<ValueType>() || keepZeros) {
                         columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue());
                         nextIndices[transition.getColumn()]++;
                     }
@@ -700,10 +714,10 @@ namespace storm {
             }
             
             storm::storage::SparseMatrix<ValueType> transposedMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), false);
-                        
+            
             return transposedMatrix;
         }
-                
+        
         template<typename ValueType>
         void SparseMatrix<ValueType>::convertToEquationSystem() {
             invertDiagonal();
@@ -713,22 +727,22 @@ namespace storm {
         template<typename ValueType>
         void SparseMatrix<ValueType>::invertDiagonal() {
             // Now iterate over all row groups and set the diagonal elements to the inverted value.
-			// If there is a row without the diagonal element, an exception is thrown.
-			ValueType one = storm::utility::one<ValueType>();
-			ValueType zero = storm::utility::zero<ValueType>();
+            // If there is a row without the diagonal element, an exception is thrown.
+            ValueType one = storm::utility::one<ValueType>();
+            ValueType zero = storm::utility::zero<ValueType>();
             bool foundDiagonalElement = false;
             for (index_type group = 0; group < this->getRowGroupCount(); ++group) {
                 for (auto& entry : this->getRowGroup(group)) {
                     if (entry.getColumn() == group) {
-						if (entry.getValue() == one) {
-							--this->nonzeroEntryCount;
-							entry.setValue(zero);
-						} else if (entry.getValue() == zero) {
-							++this->nonzeroEntryCount;
-							entry.setValue(one);
-						} else {
-							entry.setValue(one - entry.getValue());
-						}
+                        if (entry.getValue() == one) {
+                            --this->nonzeroEntryCount;
+                            entry.setValue(zero);
+                        } else if (entry.getValue() == zero) {
+                            ++this->nonzeroEntryCount;
+                            entry.setValue(one);
+                        } else {
+                            entry.setValue(one - entry.getValue());
+                        }
                         foundDiagonalElement = true;
                     }
                 }
@@ -758,7 +772,7 @@ namespace storm {
             for (index_type group = 0; group < this->getRowGroupCount(); ++group) {
                 for (auto& entry : this->getRowGroup(group)) {
                     if (entry.getColumn() == group) {
-						--this->nonzeroEntryCount;
+                        --this->nonzeroEntryCount;
                         entry.setValue(storm::utility::zero<ValueType>());
                     }
                 }
@@ -839,15 +853,28 @@ namespace storm {
             typename std::vector<ValueType>::iterator resultIterator = result.begin();
             typename std::vector<ValueType>::iterator resultIteratorEnd = result.end();
             
-            for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) {
-                *resultIterator = storm::utility::zero<ValueType>();
-                
-                for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
-                    *resultIterator += it->getValue() * vector[it->getColumn()];
+            // If the vector to multiply with and the target vector are actually the same, we need an auxiliary variable
+            // to store the intermediate result.
+            if (&vector == &result) {
+                for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) {
+                    ValueType tmpValue = storm::utility::zero<ValueType>();
+                    
+                    for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
+                        tmpValue += it->getValue() * vector[it->getColumn()];
+                    }
+                    *resultIterator = tmpValue;
+                }
+            } else {
+                for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) {
+                    *resultIterator = storm::utility::zero<ValueType>();
+                    
+                    for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
+                        *resultIterator += it->getValue() * vector[it->getColumn()];
+                    }
                 }
             }
         }
-
+        
 #ifdef STORM_HAVE_INTELTBB
         template<typename ValueType>
         void SparseMatrix<ValueType>::multiplyWithVectorParallel(std::vector<ValueType> const& vector, std::vector<ValueType>& result) const {
@@ -872,7 +899,37 @@ namespace storm {
                               });
         }
 #endif
-
+        
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::performSuccessiveOverRelaxationStep(ValueType omega, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+            const_iterator it = this->begin();
+            const_iterator ite;
+            std::vector<index_type>::const_iterator rowIterator = rowIndications.begin();
+            typename std::vector<ValueType>::const_iterator bIt = b.begin();
+            typename std::vector<ValueType>::const_iterator bIte = b.end();
+            typename std::vector<ValueType>::iterator resultIterator = x.begin();
+            typename std::vector<ValueType>::iterator resultIteratorEnd = x.end();
+            
+            // If the vector to multiply with and the target vector are actually the same, we need an auxiliary variable
+            // to store the intermediate result.
+            index_type currentRow = 0;
+            for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++bIt) {
+                ValueType tmpValue = storm::utility::zero<ValueType>();
+                ValueType diagonalElement = storm::utility::zero<ValueType>();
+                
+                for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
+                    if (it->getColumn() != currentRow) {
+                        tmpValue += it->getValue() * x[it->getColumn()];
+                    } else {
+                        diagonalElement += it->getValue();
+                    }
+                }
+                
+                *resultIterator = ((storm::utility::one<ValueType>() - omega) * *resultIterator) + (omega / diagonalElement) * (*bIt - tmpValue);
+                ++currentRow;
+            }
+        }
+        
         template<typename ValueType>
         void SparseMatrix<ValueType>::multiplyVectorWithMatrix(std::vector<value_type> const& vector, std::vector<value_type>& result) const {
             const_iterator it = this->begin();
@@ -982,7 +1039,7 @@ namespace storm {
             if (this->getRowCount() != matrix.getRowCount()) return false;
             if (this->getColumnCount() != matrix.getColumnCount()) return false;
             if (this->getRowGroupIndices() != matrix.getRowGroupIndices()) return false;
-
+            
             // Check the subset property for all rows individually.
             for (index_type row = 0; row < this->getRowCount(); ++row) {
                 for (const_iterator it1 = this->begin(row), ite1 = this->end(row), it2 = matrix.begin(row), ite2 = matrix.end(row); it1 != ite1; ++it1) {
@@ -1054,21 +1111,21 @@ namespace storm {
         }
         
         // Explicitly instantiate the entry, builder and the matrix.
-		//double
-		template class MatrixEntry<typename SparseMatrix<double>::index_type, double>;
-		template std::ostream& operator<<(std::ostream& out, MatrixEntry<uint_fast64_t, double> const& entry);
-		template class SparseMatrixBuilder<double>;
-		template class SparseMatrix<double>;
-		template std::ostream& operator<<(std::ostream& out, SparseMatrix<double> const& matrix);
-		//float
-		template class MatrixEntry<typename SparseMatrix<float>::index_type, float>;
-		template std::ostream& operator<<(std::ostream& out, MatrixEntry<uint_fast64_t, float> const& entry);
-		template class SparseMatrixBuilder<float>;
-		template class SparseMatrix<float>;
-		template std::ostream& operator<<(std::ostream& out, SparseMatrix<float> const& matrix);
-		//int
-		template class MatrixEntry<typename SparseMatrix<int>::index_type, int>;
-		template std::ostream& operator<<(std::ostream& out, MatrixEntry<uint_fast64_t, int> const& entry);
+        //double
+        template class MatrixEntry<typename SparseMatrix<double>::index_type, double>;
+        template std::ostream& operator<<(std::ostream& out, MatrixEntry<uint_fast64_t, double> const& entry);
+        template class SparseMatrixBuilder<double>;
+        template class SparseMatrix<double>;
+        template std::ostream& operator<<(std::ostream& out, SparseMatrix<double> const& matrix);
+        //float
+        template class MatrixEntry<typename SparseMatrix<float>::index_type, float>;
+        template std::ostream& operator<<(std::ostream& out, MatrixEntry<uint_fast64_t, float> const& entry);
+        template class SparseMatrixBuilder<float>;
+        template class SparseMatrix<float>;
+        template std::ostream& operator<<(std::ostream& out, SparseMatrix<float> const& matrix);
+        //int
+        template class MatrixEntry<typename SparseMatrix<int>::index_type, int>;
+        template std::ostream& operator<<(std::ostream& out, MatrixEntry<uint_fast64_t, int> const& entry);
         template class SparseMatrixBuilder<int>;
         template class SparseMatrix<int>;
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<int> const& matrix);

src/storage/SparseMatrix.h

@@ -692,6 +692,15 @@ namespace storm {
              */
             void multiplyVectorWithMatrix(std::vector<value_type> const& vector, std::vector<value_type>& result) const;
             
+            /*!
+             * Performs one step of the successive over-relaxation technique.
+             *
+             * @param omega The Omega parameter for SOR.
+             * @param x The current solution vector. The result will be written to the very same vector.
+             * @param b The 'right-hand side' of the problem.
+             */
+            void performSuccessiveOverRelaxationStep(ValueType omega, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            
             /*!
              * Multiplies the matrix with the given vector in a sequential way and writes the result to the given result
              * vector.

src/utility/solver.cpp

@@ -41,9 +41,28 @@ namespace storm {
                 return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
             }
             
+            template<typename ValueType>
+            NativeLinearEquationSolverFactory<ValueType>::NativeLinearEquationSolverFactory() {
+                switch (storm::settings::nativeEquationSolverSettings().getLinearEquationSystemMethod()) {
+                    case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Jacobi:
+                    this->method = storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod::Jacobi;
+                    break;
+                    case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel:
+                    this->method = storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod::GaussSeidel;
+                    case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::SOR:
+                    this->method = storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod::SOR;
+                }
+                omega = storm::settings::nativeEquationSolverSettings().getOmega();
+            }
+            
+            template<typename ValueType>
+            NativeLinearEquationSolverFactory<ValueType>::NativeLinearEquationSolverFactory(typename storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod method, ValueType omega) : method(method), omega(omega) {
+                // Intentionally left empty.
+            }
+            
             template<typename ValueType>
             std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
-                return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix));
+                return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix, method));
             }
             
             template<typename ValueType>

src/utility/solver.h

@@ -5,10 +5,12 @@
 
 #include "src/solver/SymbolicLinearEquationSolver.h"
 #include "src/solver/LinearEquationSolver.h"
+#include "src/solver/NativeLinearEquationSolver.h"
 #include "src/solver/MinMaxLinearEquationSolver.h"
 #include "src/solver/LpSolver.h"
 
 #include "src/storage/dd/DdType.h"
+#include "src/settings/modules/NativeEquationSolverSettings.h"
 
 #include "src/exceptions/InvalidSettingsException.h"
 
@@ -42,7 +44,14 @@ namespace storm {
             template<typename ValueType>
             class NativeLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
             public:
+                NativeLinearEquationSolverFactory();
+                NativeLinearEquationSolverFactory(typename storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod method, ValueType omega);
+                
                 virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+                
+            private:
+                typename storm::solver::NativeLinearEquationSolver<ValueType>::SolutionMethod method;
+                ValueType omega;
             };
             
             template<typename ValueType>

test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp

@@ -1,6 +1,7 @@
 #include "gtest/gtest.h"
 #include "storm-config.h"
 
+#include "src/settings/SettingMemento.h"
 #include "src/logic/Formulas.h"
 #include "src/utility/solver.h"
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
@@ -144,7 +145,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
 		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
 
-		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9)));
 
 		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
 		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
@@ -169,7 +170,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
 		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
 
-		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9)));
 
 		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
 		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
@@ -194,7 +195,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
 		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
 
-		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9)));
 
 		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
 		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
@@ -255,7 +256,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA) {
 
 		mdp.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
 
-		storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolver<double>::SolutionMethod::SOR, 0.9)));
 
 		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
 		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);