diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
index 1a1aff62f..25328e0a5 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
@@ -436,6 +436,22 @@ namespace storm {
             return result;
         }
         
+        template<class ValueType>
+        storm::storage::SparseMatrix<ValueType> SparseCtmcCslModelChecker<ValueType>::computeGeneratorMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates) {
+            storm::storage::SparseMatrix<ValueType> generatorMatrix(rateMatrix, true);
+            
+            // Place the negative exit rate on the diagonal.
+            for (uint_fast64_t row = 0; row < generatorMatrix.getRowCount(); ++row) {
+                for (auto& entry : generatorMatrix.getRow(row)) {
+                    if (entry.getColumn() == row) {
+                        entry.setValue(-exitRates[row]);
+                    }
+                }
+            }
+            
+            return generatorMatrix;
+        }
+        
         template<class ValueType>
         std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
@@ -468,9 +484,15 @@ namespace storm {
             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->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());
+            
+            std::vector<ValueType> result = SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(this->getModel(), uniformizedMatrix, subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory);
+            
             
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(this->getModel(), uniformizedMatrix, subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory)));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>());
         }
         
         // Explicitly instantiate the model checker.
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
index 529c23128..35edab8cf 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
@@ -74,8 +74,18 @@ namespace storm {
              *
              * @param rateMatrix The rate matrix.
              * @param exitRates The exit rate vector.
+             * @return The †ransition matrix of the embedded DTMC.
              */
             static storm::storage::SparseMatrix<ValueType> computeProbabilityMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates);
+
+            /*!
+             * Converts the given rate-matrix into the generator matrix.
+             *
+             * @param rateMatrix The rate matrix.
+             * @param exitRates The exit rate vector.
+             * @return The generator matrix.
+             */
+            static storm::storage::SparseMatrix<ValueType> computeGeneratorMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates);
             
             // An object that is used for solving linear equations and performing matrix-vector multiplication.
             std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index e60a6403e..7237fe765 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -325,25 +325,39 @@ namespace storm {
 
 			// 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
+            // 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)) {
@@ -353,19 +367,29 @@ namespace storm {
 				}
 			}
             
-			// Now transpose the matrix. This internally removes all explicit zeros from the matrix that were.
-            // introduced when subtracting the identity matrix.
-			bsccEquationSystem = bsccEquationSystem.transpose();
+            std::cout << bsccEquationSystem << std::endl;
+            
+            std::cout << statesInBsccsWithoutFirst << " // " << statesInBsccs << std::endl;
+            bsccEquationSystem = bsccEquationSystem.getSubmatrix(false, statesInBsccsWithoutFirst, statesInBsccs, false);
+            
+            std::cout << bsccEquationSystem << std::endl;
 
-            // Add a row to the matrix that expresses that the sum over all entries needs to be one.
+            // 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();
-            for (uint_fast64_t i = 0; i <= row; ++i) {
-                builder.addNextValue(row + 1, i, 1);
+            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;
             }
-            builder.addNextValue(row + 1, row + 1, 0);
             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);
@@ -374,6 +398,14 @@ namespace storm {
 				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.
@@ -426,12 +458,12 @@ namespace storm {
             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
+                    // 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
+                    // 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;
                 }
diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp
index 4fb0a7ed8..cbf2eafc8 100644
--- a/src/solver/GmmxxLinearEquationSolver.cpp
+++ b/src/solver/GmmxxLinearEquationSolver.cpp
@@ -149,6 +149,12 @@ 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));
diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp
index 4974720a8..2db8b576b 100644
--- a/src/storage/SparseMatrix.cpp
+++ b/src/storage/SparseMatrix.cpp
@@ -253,6 +253,13 @@ namespace storm {
             // Intentionally left empty.
         }
         
+        template<typename ValueType>
+        SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<value_type> const& other, bool insertDiagonalElements) {
+            storm::storage::BitVector rowConstraint(other.getRowCount(), true);
+            storm::storage::BitVector columnConstraint(other.getColumnCount(), true);
+            *this = other.getSubmatrix(false, rowConstraint, columnConstraint, insertDiagonalElements);
+        }
+        
         template<typename ValueType>
         SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType>&& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(std::move(other.columnsAndValues)), rowIndications(std::move(other.rowIndications)), nontrivialRowGrouping(other.nontrivialRowGrouping), rowGroupIndices(std::move(other.rowGroupIndices)) {
             // Now update the source matrix
@@ -506,9 +513,35 @@ namespace storm {
                 
         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 {
-            // First, we need to determine the number of entries and the number of rows of the submatrix.
+            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);
+            
+            // 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);
+                    ++lastIndex;
+                }
+                ++currentNumberOfSetBits;
+            }
+            
+            // Then, we need to determine the number of entries and the number of rows of the submatrix.
             index_type subEntries = 0;
             index_type subRows = 0;
+            index_type rowGroupCount = 0;
             for (auto index : rowGroupConstraint) {
                 subRows += rowGroupIndices[index + 1] - rowGroupIndices[index];
                 for (index_type i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) {
@@ -525,39 +558,18 @@ namespace storm {
                     }
                     
                     // If requested, we need to reserve one entry more for inserting the diagonal zero entry.
-                    if (insertDiagonalEntries && !foundDiagonalElement) {
+                    if (insertDiagonalEntries && !foundDiagonalElement && rowGroupCount < submatrixColumnCount) {
                         ++subEntries;
                     }
                 }
+                ++rowGroupCount;
             }
             
             // Create and initialize resulting matrix.
-            SparseMatrixBuilder<ValueType> matrixBuilder(subRows, columnConstraint.getNumberOfSetBits(), subEntries, true, this->hasNontrivialRowGrouping());
-            
-            // Create 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);
-            
-            // 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 : columnBitCountConstraint) {
-                while (lastIndex <= index) {
-                    bitsSetBeforeIndex.push_back(currentNumberOfSetBits);
-                    ++lastIndex;
-                }
-                ++currentNumberOfSetBits;
-            }
+            SparseMatrixBuilder<ValueType> matrixBuilder(subRows, submatrixColumnCount, subEntries, true, this->hasNontrivialRowGrouping());
             
             // Copy over selected entries.
+            rowGroupCount = 0;
             index_type rowCount = 0;
             for (auto index : rowGroupConstraint) {
                 if (this->hasNontrivialRowGrouping()) {
@@ -571,18 +583,18 @@ namespace storm {
                             if (index == it->getColumn()) {
                                 insertedDiagonalElement = true;
                             } else if (insertDiagonalEntries && !insertedDiagonalElement && it->getColumn() > index) {
-                                matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero<ValueType>());
+                                matrixBuilder.addNextValue(rowCount, rowCount, storm::utility::zero<ValueType>());
                                 insertedDiagonalElement = true;
                             }
                             matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[it->getColumn()], it->getValue());
                         }
                     }
-                    if (insertDiagonalEntries && !insertedDiagonalElement) {
-                        matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::zero<ValueType>());
+                    if (insertDiagonalEntries && !insertedDiagonalElement && rowGroupCount < submatrixColumnCount) {
+                        matrixBuilder.addNextValue(rowGroupCount, rowGroupCount, storm::utility::zero<ValueType>());
                     }
-                    
                     ++rowCount;
                 }
+                ++rowGroupCount;
             }
             
             return matrixBuilder.build();
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index 444c2e2f3..33863da48 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -411,6 +411,15 @@ namespace storm {
              */
             SparseMatrix(SparseMatrix<value_type> const& other);
             
+            /*!
+             * Constructs a sparse matrix by performing a deep-copy of the given matrix.
+             *
+             * @param other The matrix from which to copy the content.
+             * @param insertDiagonalElements If set to true, the copy will have all diagonal elements. If they did not
+             * exist in the original matrix, they are inserted and set to value zero.
+             */
+            SparseMatrix(SparseMatrix<value_type> const& other, bool insertDiagonalElements);
+            
             /*!
              * Constructs a sparse matrix by moving the contents of the given matrix to the newly created one.
              *
diff --git a/src/utility/vector.h b/src/utility/vector.h
index b5567ae5b..634d2bc5a 100644
--- a/src/utility/vector.h
+++ b/src/utility/vector.h
@@ -376,7 +376,9 @@ namespace storm {
 					if (val2 == 0) {
 						return (std::abs(val1) <= precision);
 					}
-                    if (std::abs((val1 - val2)/val2) > precision) return false;
+                    if (std::abs((val1 - val2)/val2) > precision) {
+                        return false;
+                    }
                 } else {
                     if (std::abs(val1 - val2) > precision) return false;
                 }
diff --git a/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp
index 08fb5a29c..9bbedfa7b 100644
--- a/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp
@@ -127,3 +127,148 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, SynchronousLeader) {
 
 	EXPECT_NEAR(1.044879046, quantitativeResult3[0], storm::settings::gmmxxEquationSolverSettings().getPrecision());
 }
+
+TEST(GmmxxDtmcPrctlModelCheckerTest, LRASingleBscc) {
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc;
+    
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 2);
+        matrixBuilder.addNextValue(0, 1, 1.);
+        matrixBuilder.addNextValue(1, 0, 1.);
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(2);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 1);
+        
+        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>()));
+        
+        auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+        auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+    }
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 4);
+        matrixBuilder.addNextValue(0, 0, .5);
+        matrixBuilder.addNextValue(0, 1, .5);
+        matrixBuilder.addNextValue(1, 0, .5);
+        matrixBuilder.addNextValue(1, 1, .5);
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(2);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 1);
+        
+        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>()));
+        
+        auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+        auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+    }
+    
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(3, 3, 3);
+        matrixBuilder.addNextValue(0, 1, 1);
+        matrixBuilder.addNextValue(1, 2, 1);
+        matrixBuilder.addNextValue(2, 0, 1);
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(3);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 2);
+        
+        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::GmmxxLinearEquationSolverFactory<double>()));
+        
+        auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+        auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+    }
+}
+
+TEST(GmmxxDtmcPrctlModelCheckerTest, LRA) {
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> mdp;
+    
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(15, 15, 20, true);
+        matrixBuilder.addNextValue(0, 1, 1);
+        matrixBuilder.addNextValue(1, 4, 0.7);
+        matrixBuilder.addNextValue(1, 6, 0.3);
+        matrixBuilder.addNextValue(2, 0, 1);
+        
+        matrixBuilder.addNextValue(3, 5, 0.8);
+        matrixBuilder.addNextValue(3, 9, 0.2);
+        matrixBuilder.addNextValue(4, 3, 1);
+        matrixBuilder.addNextValue(5, 3, 1);
+        
+        matrixBuilder.addNextValue(6, 7, 1);
+        matrixBuilder.addNextValue(7, 8, 1);
+        matrixBuilder.addNextValue(8, 6, 1);
+        
+        matrixBuilder.addNextValue(9, 10, 1);
+        matrixBuilder.addNextValue(10, 9, 1);
+        matrixBuilder.addNextValue(11, 9, 1);
+        
+        matrixBuilder.addNextValue(12, 5, 0.4);
+        matrixBuilder.addNextValue(12, 8, 0.3);
+        matrixBuilder.addNextValue(12, 11, 0.3);
+        
+        matrixBuilder.addNextValue(13, 7, 0.7);
+        matrixBuilder.addNextValue(13, 12, 0.3);
+        
+        matrixBuilder.addNextValue(14, 12, 1);
+        
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(15);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 1);
+        ap.addLabelToState("a", 4);
+        ap.addLabelToState("a", 5);
+        ap.addLabelToState("a", 7);
+        ap.addLabelToState("a", 11);
+        ap.addLabelToState("a", 13);
+        ap.addLabelToState("a", 14);
+        
+        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::GmmxxLinearEquationSolverFactory<double>()));
+        
+        auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+        auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(labelFormula);
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(0.3 / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(0.0, quantitativeResult1[3], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(1. / 3., quantitativeResult1[6], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(0.0, quantitativeResult1[9], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(0.3/3., quantitativeResult1[12], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(.79 / 3., quantitativeResult1[13], storm::settings::nativeEquationSolverSettings().getPrecision());
+        EXPECT_NEAR(0.3 / 3., quantitativeResult1[14], storm::settings::nativeEquationSolverSettings().getPrecision());
+    }
+}
diff --git a/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
index 7f53f45b4..387f54ad4 100644
--- a/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
@@ -128,7 +128,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, SynchronousLeader) {
     EXPECT_NEAR(1.0448979589010925, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 }
 
-TEST(SparseDtmcPrctlModelCheckerTest, LRA_SingleBscc) {
+TEST(SparseDtmcPrctlModelCheckerTest, LRASingleBscc) {
 	storm::storage::SparseMatrixBuilder<double> matrixBuilder;
 	std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc;