diff --git a/src/adapters/ExplicitModelAdapter.h b/src/adapters/ExplicitModelAdapter.h
index 9001803ba..aebbd1fcf 100644
--- a/src/adapters/ExplicitModelAdapter.h
+++ b/src/adapters/ExplicitModelAdapter.h
@@ -542,7 +542,7 @@ namespace storm {
                         } else {
                             // If the model is nondeterministic, we add all choices individually.
                             nondeterministicChoiceIndices.push_back(currentRow);
-                            transitionMatrixBuilder.addRowGroup();
+                            transitionMatrixBuilder.newRowGroup(currentRow);
                             
                             // First, process all unlabeled choices.
                             for (auto const& choice : allUnlabeledChoices) {
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
index 4208501d0..2f03603d8 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
@@ -122,7 +122,7 @@ public:
             // In this case we have have to compute the probabilities.
 
             // We can eliminate the rows and columns from the original transition probability matrix that have probability 0.
-            storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(statesWithProbabilityGreater0);
+            storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, statesWithProbabilityGreater0, statesWithProbabilityGreater0, true);
             
             // Compute the new set of target states in the reduced system.
             storm::storage::BitVector rightStatesInReducedSystem = psiStates % statesWithProbabilityGreater0;
@@ -291,7 +291,7 @@ public:
             // In this case we have have to compute the probabilities.
             
             // We can eliminate the rows and columns from the original transition probability matrix.
-            storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates);
+            storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates, true);
             // Converting the matrix from the fixpoint notation to the form needed for the equation
             // system. That is, we go from x = A*x + b to (I-A)x = b.
             submatrix.convertToEquationSystem();
@@ -444,7 +444,7 @@ public:
         } else {
             // In this case we have to compute the reward values for the remaining states.
             // We can eliminate the rows and columns from the original transition probability matrix.
-            storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates);
+            storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates, true);
             // Converting the matrix from the fixpoint notation to the form needed for the equation
             // system. That is, we go from x = A*x + b to (I-A)x = b.
             submatrix.convertToEquationSystem();
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
index 7aa5ae581..0e319e67f 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
@@ -117,7 +117,7 @@ namespace storm {
                         // In this case we have have to compute the probabilities.
                         
                         // We can eliminate the rows and columns from the original transition probability matrix that have probability 0.
-                        storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(statesWithProbabilityGreater0, this->getModel().getNondeterministicChoiceIndices());
+                        storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, statesWithProbabilityGreater0, statesWithProbabilityGreater0, false);
                         
                         // Get the "new" nondeterministic choice indices for the submatrix.
                         std::vector<uint_fast64_t> subNondeterministicChoiceIndices = storm::utility::vector::getConstrainedOffsetVector(this->getModel().getNondeterministicChoiceIndices(), statesWithProbabilityGreater0);
@@ -126,7 +126,7 @@ namespace storm {
                         storm::storage::BitVector rightStatesInReducedSystem = psiStates % statesWithProbabilityGreater0;
                         
                         // Make all rows absorbing that satisfy the second sub-formula.
-                        submatrix.makeRowsAbsorbing(rightStatesInReducedSystem, subNondeterministicChoiceIndices);
+                        submatrix.makeRowGroupsAbsorbing(rightStatesInReducedSystem, subNondeterministicChoiceIndices);
                         
                         // Create the vector with which to multiply.
                         std::vector<Type> subresult(statesWithProbabilityGreater0.getNumberOfSetBits());
@@ -320,14 +320,14 @@ namespace storm {
 
                         // First, we can eliminate the rows and columns from the original transition probability matrix for states
                         // whose probabilities are already known.
-                        storm::storage::SparseMatrix<Type> submatrix = transitionMatrix.getSubmatrix(maybeStates, nondeterministicChoiceIndices);
+                        storm::storage::SparseMatrix<Type> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, false);
                         
                         // Get the "new" nondeterministic choice indices for the submatrix.
                         std::vector<uint_fast64_t> subNondeterministicChoiceIndices = storm::utility::vector::getConstrainedOffsetVector(nondeterministicChoiceIndices, maybeStates);
                         
                         // Prepare the right-hand side of the equation system. For entry i this corresponds to
                         // the accumulated probability of going from state i to some 'yes' state.
-                        std::vector<Type> b = transitionMatrix.getConstrainedRowSumVector(maybeStates, nondeterministicChoiceIndices, statesWithProbability1);
+                        std::vector<Type> b = transitionMatrix.getConstrainedRowGroupSumVector(maybeStates, nondeterministicChoiceIndices, statesWithProbability1);
                         
                         // Create vector for results for maybe states.
                         std::vector<Type> x(maybeStates.getNumberOfSetBits());
@@ -487,7 +487,7 @@ namespace storm {
                         
                         // We can eliminate the rows and columns from the original transition probability matrix for states
                         // whose reward values are already known.
-                        storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates, this->getModel().getNondeterministicChoiceIndices());
+                        storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates, false);
                         
                         // Get the "new" nondeterministic choice indices for the submatrix.
                         std::vector<uint_fast64_t> subNondeterministicChoiceIndices = storm::utility::vector::getConstrainedOffsetVector(this->getModel().getNondeterministicChoiceIndices(), maybeStates);
diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp
index a2a5e6b82..cb62d529c 100644
--- a/src/storage/SparseMatrix.cpp
+++ b/src/storage/SparseMatrix.cpp
@@ -71,7 +71,7 @@ namespace storm {
                 
                 if (!hasCustomRowGrouping) {
                     for (uint_fast64_t i = lastRow + 1; i <= row; ++i) {
-                        rowGroupIndices.push_back(row);
+                        rowGroupIndices.push_back(i);
                         ++currentRowGroup;
                     }
                 }
@@ -97,15 +97,15 @@ namespace storm {
         }
         
         template<typename T>
-        void SparseMatrixBuilder<T>::addRowGroup() {
+        void SparseMatrixBuilder<T>::newRowGroup(uint_fast64_t startingRow) {
             if (!hasCustomRowGrouping) {
                 throw storm::exceptions::InvalidStateException() << "Illegal call to SparseMatrix::addRowGroup: matrix was not created to have a custom row grouping.";
             }
             if (rowGroupCountSet) {
-                rowGroupIndices[currentRowGroup] = lastRow + 1;
+                rowGroupIndices[currentRowGroup] = startingRow;
                 ++currentRowGroup;
             } else {
-                rowGroupIndices.push_back(lastRow + 1);
+                rowGroupIndices.push_back(startingRow);
             }
         }
         
@@ -282,10 +282,10 @@ namespace storm {
             }
             
             bool equalityResult = true;
-            
-            equalityResult &= rowCount == other.rowCount;
-            equalityResult &= columnCount == other.columnCount;
-            equalityResult &= rowGroupIndices == other.rowGroupIndices;
+
+            equalityResult &= this->getRowCount() == other.getRowCount();
+            equalityResult &= this->getColumnCount() == other.getColumnCount();
+            equalityResult &= this->getRowGroupIndices() == other.getRowGroupIndices();
             
             // For the actual contents, we need to do a little bit more work, because we want to ignore elements that
             // are set to zero, please they may be represented implicitly in the other matrix.
@@ -333,28 +333,29 @@ namespace storm {
             return rowGroupIndices.size() - 1;
         }
         
+        template<typename T>
+        std::vector<uint_fast64_t> const& SparseMatrix<T>::getRowGroupIndices() const {
+            return rowGroupIndices;
+        }
+
         template<typename T>
         void SparseMatrix<T>::makeRowsAbsorbing(storm::storage::BitVector const& rows) {
             for (auto row : rows) {
-                makeRowAbsorbing(row, row);
+                makeRowDirac(row, row);
             }
         }
         
         template<typename T>
-        void SparseMatrix<T>::makeRowsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices) {
+        void SparseMatrix<T>::makeRowGroupsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices) {
             for (auto rowGroup : rowGroupConstraint) {
-                for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) {
-                    makeRowAbsorbing(row, rowGroup);
+                for (uint_fast64_t row = this->getRowGroupIndices()[rowGroup]; row < this->getRowGroupIndices()[rowGroup + 1]; ++row) {
+                    makeRowDirac(row, rowGroup);
                 }
             }
         }
         
         template<typename T>
-        void SparseMatrix<T>::makeRowAbsorbing(uint_fast64_t row, uint_fast64_t column) {
-            if (row > rowCount) {
-                throw storm::exceptions::OutOfRangeException() << "Illegal call to SparseMatrix::makeRowAbsorbing: access to row " << row << " is out of bounds.";
-            }
-            
+        void SparseMatrix<T>::makeRowDirac(uint_fast64_t row, uint_fast64_t column) {
             iterator columnValuePtr = this->begin(row);
             iterator columnValuePtrEnd = this->end(row);
             
@@ -397,11 +398,11 @@ namespace storm {
         }
         
         template<typename T>
-        std::vector<T> SparseMatrix<T>::getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, storm::storage::BitVector const& columnConstraint) const {
+        std::vector<T> SparseMatrix<T>::getConstrainedRowGroupSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, storm::storage::BitVector const& columnConstraint) const {
             std::vector<T> result;
             result.reserve(rowGroupConstraint.getNumberOfSetBits());
             for (auto rowGroup : rowGroupConstraint) {
-                for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) {
+                for (uint_fast64_t row = this->getRowGroupIndices()[rowGroup]; row < this->getRowGroupIndices()[rowGroup + 1]; ++row) {
                     result.push_back(getConstrainedRowSum(row, columnConstraint));
                 }
             }
@@ -409,21 +410,20 @@ namespace storm {
         }
         
         template<typename T>
-        SparseMatrix<T> SparseMatrix<T>::getSubmatrix(storm::storage::BitVector const& constraint) const {
-            // Create a fake row grouping to reduce this to a call to a more general method.
-            std::vector<uint_fast64_t> rowGroupIndices(rowCount + 1);
-            uint_fast64_t i = 0;
-            for (std::vector<uint_fast64_t>::iterator it = rowGroupIndices.begin(); it != rowGroupIndices.end(); ++it, ++i) {
-                *it = i;
+        SparseMatrix<T> SparseMatrix<T>::getSubmatrix(bool useGroups, storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint, bool insertDiagonalElements) const {
+            if (useGroups) {
+                return getSubmatrix(rowConstraint, columnConstraint, this->getRowGroupIndices(), insertDiagonalElements);
+            } else {
+                // Create a fake row grouping to reduce this to a call to a more general method.
+                std::vector<uint_fast64_t> fakeRowGroupIndices(rowCount + 1);
+                uint_fast64_t i = 0;
+                for (std::vector<uint_fast64_t>::iterator it = fakeRowGroupIndices.begin(); it != fakeRowGroupIndices.end(); ++it, ++i) {
+                    *it = i;
+                }
+                return getSubmatrix(rowConstraint, columnConstraint, fakeRowGroupIndices, insertDiagonalElements);
             }
-            return getSubmatrix(constraint, constraint, rowGroupIndices);
         }
-        
-        template<typename T>
-        SparseMatrix<T> SparseMatrix<T>::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries) const {
-            return getSubmatrix(rowGroupConstraint, rowGroupConstraint, rowGroupIndices, insertDiagonalEntries);
-        }
-        
+                
         template<typename T>
         SparseMatrix<T> SparseMatrix<T>::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries) const {
             // First, we need to determine the number of entries and the number of rows of the submatrix.
@@ -452,7 +452,7 @@ namespace storm {
             }
             
             // Create and initialize resulting matrix.
-            SparseMatrixBuilder<T> matrixBuilder(subRows, columnConstraint.getNumberOfSetBits(), subEntries);
+            SparseMatrixBuilder<T> matrixBuilder(subRows, columnConstraint.getNumberOfSetBits(), subEntries, true);
             
             // 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.
@@ -480,6 +480,7 @@ namespace storm {
             // Copy over selected entries.
             uint_fast64_t rowCount = 0;
             for (auto index : rowGroupConstraint) {
+                matrixBuilder.newRowGroup(rowCount);
                 for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) {
                     bool insertedDiagonalElement = false;
                     
@@ -506,7 +507,7 @@ namespace storm {
         }
         
         template<typename T>
-        SparseMatrix<T> SparseMatrix<T>::getSubmatrix(std::vector<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries) const {
+        SparseMatrix<T> SparseMatrix<T>::selectRowsFromRowGroups(std::vector<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries) const {
             // First, we need to count how many non-zero entries the resulting matrix will have and reserve space for
             // diagonal entries if requested.
             uint_fast64_t subEntries = 0;
@@ -612,20 +613,14 @@ namespace storm {
         
         template<typename T>
         void SparseMatrix<T>::invertDiagonal() {
-            // Check if the matrix is square, because only then it makes sense to perform this
-            // transformation.
-            if (this->getRowCount() != this->getColumnCount()) {
-                throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to be square!";
-            }
-            
-            // Now iterate over all rows and set the diagonal elements to the inverted value.
+            // 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.
             T one = storm::utility::constantOne<T>();
             bool foundDiagonalElement = false;
-            for (uint_fast64_t row = 0; row < rowCount; ++row) {
-                for (iterator it = this->begin(row), ite = this->end(row); it != ite; ++it) {
-                    if (it->first == row) {
-                        it->second = one - it->second;
+            for (uint_fast64_t group = 0; group < this->getRowGroupCount(); ++group) {
+                for (auto& entry : this->getRowGroup(group)) {
+                    if (entry.first == group) {
+                        entry.second = one - entry.second;
                         foundDiagonalElement = true;
                     }
                 }
@@ -639,16 +634,11 @@ namespace storm {
         
         template<typename T>
         void SparseMatrix<T>::negateAllNonDiagonalEntries() {
-            // Check if the matrix is square, because only then it makes sense to perform this transformation.
-            if (this->getRowCount() != this->getColumnCount()) {
-                throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrix::invertDiagonal: matrix is non-square.";
-            }
-            
-            // Iterate over all rows and negate all the elements that are not on the diagonal.
-            for (uint_fast64_t row = 0; row < rowCount; ++row) {
-                for (iterator it = this->begin(row), ite = this->end(row); it != ite; ++it) {
-                    if (it->first != row) {
-                        it->second = -it->second;
+            // Iterate over all row groups and negate all the elements that are not on the diagonal.
+            for (uint_fast64_t group = 0; group < this->getRowGroupCount(); ++group) {
+                for (auto& entry : this->getRowGroup(group)) {
+                    if (entry.first != group) {
+                        entry.second = -entry.second;
                     }
                 }
             }
@@ -656,16 +646,11 @@ namespace storm {
         
         template<typename T>
         void SparseMatrix<T>::deleteDiagonalEntries() {
-            // Check if the matrix is square, because only then it makes sense to perform this transformation.
-            if (this->getRowCount() != this->getColumnCount()) {
-                throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrix::deleteDiagonalEntries: matrix is non-square.";
-            }
-            
             // Iterate over all rows and negate all the elements that are not on the diagonal.
-            for (uint_fast64_t row = 0; row < rowCount; ++row) {
-                for (iterator it = this->begin(row), ite = this->end(row); it != ite; ++it) {
-                    if (it->first == row) {
-                        it->second = storm::utility::constantZero<T>();
+            for (uint_fast64_t group = 0; group < this->getRowGroupCount(); ++group) {
+                for (auto& entry : this->getRowGroup(group)) {
+                    if (entry.first == group) {
+                        entry.second = storm::utility::constantZero<T>();
                     }
                 }
             }
@@ -853,6 +838,8 @@ 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 (uint_fast64_t 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) {
@@ -872,33 +859,36 @@ namespace storm {
         std::ostream& operator<<(std::ostream& out, SparseMatrix<T> const& matrix) {
             // Print column numbers in header.
             out << "\t\t";
-            for (uint_fast64_t i = 0; i < matrix.columnCount; ++i) {
+            for (uint_fast64_t i = 0; i < matrix.getColumnCount(); ++i) {
                 out << i << "\t";
             }
             out << std::endl;
             
-            // Iterate over all rows.
-            for (uint_fast64_t i = 0; i < matrix.rowCount; ++i) {
-                uint_fast64_t nextIndex = matrix.rowIndications[i];
-                
-                // Print the actual row.
-                out << i << "\t(\t";
-                uint_fast64_t currentRealIndex = 0;
-                while (currentRealIndex < matrix.columnCount) {
-                    if (nextIndex < matrix.rowIndications[i + 1] && currentRealIndex == matrix.columnsAndValues[nextIndex].first) {
-                        out << matrix.columnsAndValues[nextIndex].second << "\t";
-                        ++nextIndex;
-                    } else {
-                        out << "0\t";
+            // Iterate over all row groups.
+            for (uint_fast64_t group = 0; group < matrix.getRowGroupCount(); ++group) {
+                out << "\t---- group " << group << " ----" << std::endl;
+                for (uint_fast64_t i = matrix.getRowGroupIndices()[group]; i < matrix.getRowGroupIndices()[group + 1]; ++i) {
+                    uint_fast64_t nextIndex = matrix.rowIndications[i];
+                    
+                    // Print the actual row.
+                    out << i << "\t(\t";
+                    uint_fast64_t currentRealIndex = 0;
+                    while (currentRealIndex < matrix.columnCount) {
+                        if (nextIndex < matrix.rowIndications[i + 1] && currentRealIndex == matrix.columnsAndValues[nextIndex].first) {
+                            out << matrix.columnsAndValues[nextIndex].second << "\t";
+                            ++nextIndex;
+                        } else {
+                            out << "0\t";
+                        }
+                        ++currentRealIndex;
                     }
-                    ++currentRealIndex;
+                    out << "\t)\t" << i << std::endl;
                 }
-                out << "\t)\t" << i << std::endl;
             }
             
             // Print column numbers in footer.
             out << "\t\t";
-            for (uint_fast64_t i = 0; i < matrix.columnCount; ++i) {
+            for (uint_fast64_t i = 0; i < matrix.getColumnCount(); ++i) {
                 out << i << "\t";
             }
             out << std::endl;
@@ -910,11 +900,12 @@ namespace storm {
         std::size_t SparseMatrix<T>::hash() const {
             std::size_t result = 0;
             
-            boost::hash_combine(result, rowCount);
-            boost::hash_combine(result, columnCount);
-            boost::hash_combine(result, entryCount);
+            boost::hash_combine(result, this->getRowCount());
+            boost::hash_combine(result, this->getColumnCount());
+            boost::hash_combine(result, this->getEntryCount());
             boost::hash_combine(result, boost::hash_range(columnsAndValues.begin(), columnsAndValues.end()));
             boost::hash_combine(result, boost::hash_range(rowIndications.begin(), rowIndications.end()));
+            boost::hash_combine(result, boost::hash_range(rowGroupIndices.begin(), rowGroupIndices.end()));
             
             return result;
         }
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index f066744de..83eb4e46d 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -62,10 +62,12 @@ namespace storm {
             void addNextValue(uint_fast64_t row, uint_fast64_t column, T const& value);
             
             /*!
-             * Adds a new row group to the matrix. Note that this needs to be called before any entries in the new row
+             * Starts a new row group in the matrix. Note that this needs to be called before any entries in the new row
              * group are added.
+             *
+             * @param startingRow The starting row of the new row group.
              */
-            void addRowGroup();
+            void newRowGroup(uint_fast64_t startingRow);
             
             /*
              * Finalizes the sparse matrix to indicate that initialization process has been completed and the matrix
@@ -344,6 +346,13 @@ namespace storm {
              */
             uint_fast64_t getRowGroupCount() const;
             
+            /*!
+             * Returns the grouping of rows of this matrix.
+             *
+             * @return The grouping of rows of this matrix.
+             */
+            std::vector<uint_fast64_t> const& getRowGroupIndices() const;
+            
             /*!
              * This function makes the given rows absorbing.
              *
@@ -357,16 +366,16 @@ namespace storm {
              * @param rowGroupConstraint A bit vector indicating which row groups to make absorbing.
              * @param rowGroupIndices A vector indicating which rows belong to a given row group.
              */
-            void makeRowsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices);
+            void makeRowGroupsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices);
             
             /*!
-             * This function makes the given row absorbing. This means that all entries will be set to 0 except the one
+             * This function makes the given row Dirac. This means that all entries will be set to 0 except the one
              * at the specified column, which is set to 1 instead.
              *
-             * @param row The row to be made absorbing.
+             * @param row The row to be made Dirac.
              * @param column The index of the column whose value is to be set to 1.
              */
-            void makeRowAbsorbing(const uint_fast64_t row, const uint_fast64_t column);
+            void makeRowDirac(const uint_fast64_t row, const uint_fast64_t column);
             
             /*
              * Sums the entries in the given row and columns.
@@ -398,54 +407,32 @@ namespace storm {
              * @return A vector whose entries represent the sums of selected columns for all rows in selected row
              * groups.
              */
-            std::vector<T> getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, storm::storage::BitVector const& columnConstraint) const;
+            std::vector<T> getConstrainedRowGroupSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, storm::storage::BitVector const& columnConstraint) const;
             
             /*!
              * Creates a submatrix of the current matrix by dropping all rows and columns whose bits are not
              * set to one in the given bit vector.
              *
-             * @param constraint A bit vector indicating which rows and columns to keep.
-             * @return A matrix corresponding to a submatrix of the current matrix in which only rows and columns given
-             * by the constraint are kept and all others are dropped.
-             */
-            SparseMatrix getSubmatrix(storm::storage::BitVector const& constraint) const;
-            
-            /*!
-             * Creates a submatrix of the current matrix by keeping only certain row groups and columns.
-             *
-             * @param rowGroupConstraint A bit vector indicating which row groups and columns to keep.
-             * @param rowGroupIndices A vector indicating which rows belong to a given row group.
-             * @param insertDiagonalEntries If set to true, the resulting matrix will have zero entries in column i for
-             * each row in row group i. This can then be used for inserting other values later.
-             * @return A matrix corresponding to a submatrix of the current matrix in which only row groups and columns
-             * given by the row group constraint are kept and all others are dropped.
-             */
-            SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries = false) const;
-            
-            /*!
-             * Creates a submatrix of the current matrix by keeping only row groups and columns in the given row group
-             * and column constraint, respectively.
-             *
-             * @param rowGroupConstraint A bit vector indicating which row groups to keep.
+             * @param useGroups If set to true, the constraint for the rows is interpreted as selecting whole row groups.
+             * @param constraint A bit vector indicating which rows to keep.
              * @param columnConstraint A bit vector indicating which columns to keep.
-             * @param rowGroupIndices A vector indicating which rows belong to a given row group.
              * @param insertDiagonalEntries If set to true, the resulting matrix will have zero entries in column i for
-             * each row in row group i. This can then be used for inserting other values later.
-             * @return A matrix corresponding to a submatrix of the current matrix in which only row groups and columns
-             * given by the row group constraint are kept and all others are dropped.
+             * each row i, if there is no value yet. This can then be used for inserting other values later.
+             * @return A matrix corresponding to a submatrix of the current matrix in which only rows and columns given
+             * by the constraints are kept and all others are dropped.
              */
-            SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries = false) const;
+            SparseMatrix getSubmatrix(bool useGroups, storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint, bool insertDiagonalEntries = false) const;
             
             /*!
-             * Creates a submatrix of the current matrix by selecting one row out of each row group.
+             * Selects exactly one row from each row group of this matrix and returns the resulting matrix.
              *
-             * @param rowGroupdToRowIndexMapping A mapping from each row group index to a selected row in this group.
+             * @param rowGroupToRowIndexMapping A mapping from each row group index to a selected row in this group.
              * @param rowGroupIndices A vector indicating which rows belong to a given row group.
              * @param insertDiagonalEntries If set to true, the resulting matrix will have zero entries in column i for
              * each row in row group i. This can then be used for inserting other values later.
              * @return A submatrix of the current matrix by selecting one row out of each row group.
              */
-            SparseMatrix getSubmatrix(std::vector<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries = true) const;
+            SparseMatrix selectRowsFromRowGroups(std::vector<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries = true) const;
             
             /*!
              * Transposes the matrix.
@@ -505,6 +492,40 @@ namespace storm {
              */
             void multiplyWithVector(std::vector<T> const& vector, std::vector<T>& result) const;
             
+            /*!
+             * Computes the sum of the entries in a given row.
+             *
+             * @param row The row that is to be summed.
+             * @return The sum of the selected row.
+             */
+            T getRowSum(uint_fast64_t row) const;
+            
+            /*!
+             * Checks if the current matrix is a submatrix of the given matrix, where a matrix A is called a submatrix
+             * of B if B has no entries in position where A has none. Additionally, the matrices must be of equal size.
+             *
+             * @param matrix The matrix that possibly is a supermatrix of the current matrix.
+             * @return True iff the current matrix is a submatrix of the given matrix.
+             */
+            bool isSubmatrixOf(SparseMatrix<T> const& matrix) const;
+            
+            template<typename TPrime>
+            friend std::ostream& operator<<(std::ostream& out, SparseMatrix<TPrime> const& matrix);
+            
+            /*!
+             * Returns the size of the matrix in memory measured in bytes.
+             *
+             * @return The size of the matrix in memory measured in bytes.
+             */
+            uint_fast64_t getSizeInMemory() const;
+            
+            /*!
+             * Calculates a hash value over all values contained in the matrix.
+             *
+             * @return size_t A hash value for this matrix.
+             */
+            std::size_t hash() const;
+            
             /*!
              * Returns an object representing the consecutive rows given by the parameters.
              *
@@ -600,42 +621,22 @@ namespace storm {
              * @return An iterator that points past the end of the last row of the matrix.
              */
             iterator end();
-
-            /*!
-             * Computes the sum of the entries in a given row.
-             *
-             * @param row The row that is to be summed.
-             * @return The sum of the selected row.
-             */
-            T getRowSum(uint_fast64_t row) const;
-            
-            /*!
-             * Checks if the current matrix is a submatrix of the given matrix, where a matrix A is called a submatrix
-             * of B if B has no entries in position where A has none. Additionally, the matrices must be of equal size.
-             *
-             * @param matrix The matrix that possibly is a supermatrix of the current matrix.
-             * @return True iff the current matrix is a submatrix of the given matrix.
-             */
-            bool isSubmatrixOf(SparseMatrix<T> const& matrix) const;
-            
-            template<typename TPrime>
-            friend std::ostream& operator<<(std::ostream& out, SparseMatrix<TPrime> const& matrix);
-            
-            /*!
-             * Returns the size of the matrix in memory measured in bytes.
-             *
-             * @return The size of the matrix in memory measured in bytes.
-             */
-            uint_fast64_t getSizeInMemory() const;
             
+        private:
             /*!
-             * Calculates a hash value over all values contained in the matrix.
+             * Creates a submatrix of the current matrix by keeping only row groups and columns in the given row group
+             * and column constraint, respectively.
              *
-             * @return size_t A hash value for this matrix.
+             * @param rowGroupConstraint A bit vector indicating which row groups to keep.
+             * @param columnConstraint A bit vector indicating which columns to keep.
+             * @param rowGroupIndices A vector indicating which rows belong to a given row group.
+             * @param insertDiagonalEntries If set to true, the resulting matrix will have zero entries in column i for
+             * each row in row group i. This can then be used for inserting other values later.
+             * @return A matrix corresponding to a submatrix of the current matrix in which only row groups and columns
+             * given by the row group constraint are kept and all others are dropped.
              */
-            std::size_t hash() const;
+            SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries = false) const;
             
-        private:
             // The number of rows of the matrix.
             uint_fast64_t rowCount;
             
diff --git a/test/functional/storage/SparseMatrixTest.cpp b/test/functional/storage/SparseMatrixTest.cpp
index 88aa2cf3c..b2a53c566 100644
--- a/test/functional/storage/SparseMatrixTest.cpp
+++ b/test/functional/storage/SparseMatrixTest.cpp
@@ -221,13 +221,16 @@ TEST(SparseMatrix, MakeAbsorbing) {
 }
 
 TEST(SparseMatrix, MakeRowGroupAbsorbing) {
-    storm::storage::SparseMatrixBuilder<double> matrixBuilder(5, 4, 9);
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder(5, 4, 9, true);
+    ASSERT_NO_THROW(matrixBuilder.newRowGroup(0));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(0, 1, 1.0));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(0, 2, 1.2));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(1, 0, 0.5));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(1, 1, 0.7));
+    ASSERT_NO_THROW(matrixBuilder.newRowGroup(2));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(2, 0, 0.5));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(3, 2, 1.1));
+    ASSERT_NO_THROW(matrixBuilder.newRowGroup(4));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(4, 0, 0.1));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(4, 1, 0.2));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(4, 3, 0.3));
@@ -239,15 +242,18 @@ TEST(SparseMatrix, MakeRowGroupAbsorbing) {
     storm::storage::BitVector absorbingRowGroups(3);
     absorbingRowGroups.set(1);
     
-    ASSERT_NO_THROW(matrix.makeRowsAbsorbing(absorbingRowGroups, rowGroupIndices));
+    ASSERT_NO_THROW(matrix.makeRowGroupsAbsorbing(absorbingRowGroups, rowGroupIndices));
 
-    storm::storage::SparseMatrixBuilder<double> matrixBuilder2;
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder2(0, 0, 0, true);
+    ASSERT_NO_THROW(matrixBuilder2.newRowGroup(0));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(0, 1, 1.0));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(0, 2, 1.2));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(1, 0, 0.5));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(1, 1, 0.7));
+    ASSERT_NO_THROW(matrixBuilder2.newRowGroup(2));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(2, 1, 1));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(3, 1, 1));
+    ASSERT_NO_THROW(matrixBuilder2.newRowGroup(4));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(4, 0, 0.1));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(4, 1, 0.2));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(4, 3, 0.3));
@@ -279,14 +285,17 @@ TEST(SparseMatrix, ConstrainedRowSumVector) {
     std::vector<double> constrainedRowSum = matrix.getConstrainedRowSumVector(storm::storage::BitVector(5, true), columnConstraint);
     ASSERT_TRUE(constrainedRowSum == std::vector<double>({1.0, 0.7, 0, 0, 0.5}));
     
-    storm::storage::SparseMatrixBuilder<double> matrixBuilder2(5, 4, 9);
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder2(5, 4, 9, true);
+    ASSERT_NO_THROW(matrixBuilder2.newRowGroup(0));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(0, 1, 1.0));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(0, 2, 1.2));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(1, 0, 0.5));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(1, 1, 0.7));
+    ASSERT_NO_THROW(matrixBuilder2.newRowGroup(2));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(2, 0, 0.5));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(3, 2, 1.1));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(3, 3, 1.2));
+    ASSERT_NO_THROW(matrixBuilder2.newRowGroup(4));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(4, 0, 0.1));
     ASSERT_NO_THROW(matrixBuilder2.addNextValue(4, 3, 0.3));
     storm::storage::SparseMatrix<double> matrix2;
@@ -301,19 +310,23 @@ TEST(SparseMatrix, ConstrainedRowSumVector) {
     columnConstraint2.set(2);
     columnConstraint2.set(3);
     
-    ASSERT_NO_THROW(std::vector<double> constrainedRowSum2 = matrix2.getConstrainedRowSumVector(rowGroupConstraint, rowGroupIndices, columnConstraint2));
-    std::vector<double> constrainedRowSum2 = matrix2.getConstrainedRowSumVector(rowGroupConstraint, rowGroupIndices, columnConstraint2);
+    ASSERT_NO_THROW(std::vector<double> constrainedRowSum2 = matrix2.getConstrainedRowGroupSumVector(rowGroupConstraint, rowGroupIndices, columnConstraint2));
+    std::vector<double> constrainedRowSum2 = matrix2.getConstrainedRowGroupSumVector(rowGroupConstraint, rowGroupIndices, columnConstraint2);
     ASSERT_TRUE(constrainedRowSum2 == std::vector<double>({0, 2.3}));
 }
 
 TEST(SparseMatrix, Submatrix) {
-    storm::storage::SparseMatrixBuilder<double> matrixBuilder(5, 4, 9);
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder(5, 4, 9, true);
+    ASSERT_NO_THROW(matrixBuilder.newRowGroup(0));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(0, 1, 1.0));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(0, 2, 1.2));
+    ASSERT_NO_THROW(matrixBuilder.newRowGroup(1));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(1, 0, 0.5));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(1, 1, 0.7));
+    ASSERT_NO_THROW(matrixBuilder.newRowGroup(2));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(2, 0, 0.5));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(3, 2, 1.1));
+    ASSERT_NO_THROW(matrixBuilder.newRowGroup(4));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(4, 0, 0.1));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(4, 1, 0.2));
     ASSERT_NO_THROW(matrixBuilder.addNextValue(4, 3, 0.3));
@@ -330,11 +343,13 @@ TEST(SparseMatrix, Submatrix) {
     columnConstraint.set(0);
     columnConstraint.set(3);
     
-    ASSERT_NO_THROW(storm::storage::SparseMatrix<double> matrix2 = matrix.getSubmatrix(rowGroupConstraint, columnConstraint, rowGroupIndices, false));
-    storm::storage::SparseMatrix<double> matrix2 = matrix.getSubmatrix(rowGroupConstraint, columnConstraint, rowGroupIndices, false);
+    ASSERT_NO_THROW(storm::storage::SparseMatrix<double> matrix2 = matrix.getSubmatrix(true, rowGroupConstraint, columnConstraint, false));
+    storm::storage::SparseMatrix<double> matrix2 = matrix.getSubmatrix(true, rowGroupConstraint, columnConstraint, false);
     
-    storm::storage::SparseMatrixBuilder<double> matrixBuilder3(3, 2, 3);
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder3(3, 2, 3, true);
+    ASSERT_NO_THROW(matrixBuilder3.newRowGroup(0));
     ASSERT_NO_THROW(matrixBuilder3.addNextValue(0, 0, 0.5));
+    ASSERT_NO_THROW(matrixBuilder3.newRowGroup(2));
     ASSERT_NO_THROW(matrixBuilder3.addNextValue(2, 0, 0.1));
     ASSERT_NO_THROW(matrixBuilder3.addNextValue(2, 1, 0.3));
     storm::storage::SparseMatrix<double> matrix3;
@@ -344,8 +359,8 @@ TEST(SparseMatrix, Submatrix) {
     
     std::vector<uint_fast64_t> rowGroupToIndexMapping = {0, 0, 1, 0};
 
-    ASSERT_NO_THROW(storm::storage::SparseMatrix<double> matrix4 = matrix.getSubmatrix(rowGroupToIndexMapping, rowGroupIndices));
-    storm::storage::SparseMatrix<double> matrix4 = matrix.getSubmatrix(rowGroupToIndexMapping, rowGroupIndices);
+    ASSERT_NO_THROW(storm::storage::SparseMatrix<double> matrix4 = matrix.selectRowsFromRowGroups(rowGroupToIndexMapping, rowGroupIndices));
+    storm::storage::SparseMatrix<double> matrix4 = matrix.selectRowsFromRowGroups(rowGroupToIndexMapping, rowGroupIndices);
     
     storm::storage::SparseMatrixBuilder<double> matrixBuilder5(4, 4, 8);
     ASSERT_NO_THROW(matrixBuilder5.addNextValue(0, 1, 1.0));