diff --git a/src/builder/ExplicitPrismModelBuilder.cpp b/src/builder/ExplicitPrismModelBuilder.cpp
index 9a0408a3c..1aeab458d 100644
--- a/src/builder/ExplicitPrismModelBuilder.cpp
+++ b/src/builder/ExplicitPrismModelBuilder.cpp
@@ -24,6 +24,7 @@
 
 namespace storm {
     namespace builder {
+       
         /*!
          * A structure that is used to keep track of a reward model currently being built.
          */
@@ -82,17 +83,17 @@ namespace storm {
         }
         
         template <typename ValueType, typename RewardModelType, typename StateType>
-        ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::Options::Options() : buildCommandLabels(false), buildAllRewardModels(true), buildStateInformation(false), rewardModelsToBuild(), constantDefinitions(), buildAllLabels(true), labelsToBuild(), expressionLabels(), terminalStates(), negatedTerminalStates() {
+        ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::Options::Options() : explorationOrder(storm::settings::generalSettings().getExplorationOrder()), buildCommandLabels(false), buildAllRewardModels(true), buildStateInformation(false), rewardModelsToBuild(), constantDefinitions(), buildAllLabels(true), labelsToBuild(), expressionLabels(), terminalStates(), negatedTerminalStates() {
             // Intentionally left empty.
         }
         
         template <typename ValueType, typename RewardModelType, typename StateType>
-        ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::Options::Options(storm::logic::Formula const& formula) : buildCommandLabels(false), buildAllRewardModels(false), buildStateInformation(false), rewardModelsToBuild(), constantDefinitions(), buildAllLabels(false), labelsToBuild(std::set<std::string>()), expressionLabels(std::vector<storm::expressions::Expression>()), terminalStates(), negatedTerminalStates() {
+        ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::Options::Options(storm::logic::Formula const& formula) : explorationOrder(storm::settings::generalSettings().getExplorationOrder()), buildCommandLabels(false), buildAllRewardModels(false), buildStateInformation(false), rewardModelsToBuild(), constantDefinitions(), buildAllLabels(false), labelsToBuild(std::set<std::string>()), expressionLabels(std::vector<storm::expressions::Expression>()), terminalStates(), negatedTerminalStates() {
             this->preserveFormula(formula);
         }
         
         template <typename ValueType, typename RewardModelType, typename StateType>
-        ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::Options::Options(std::vector<std::shared_ptr<const storm::logic::Formula>> const& formulas) : buildCommandLabels(false), buildAllRewardModels(false), buildStateInformation(false), rewardModelsToBuild(), constantDefinitions(), buildAllLabels(false), labelsToBuild(), expressionLabels(), terminalStates(), negatedTerminalStates() {
+        ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::Options::Options(std::vector<std::shared_ptr<const storm::logic::Formula>> const& formulas) : explorationOrder(storm::settings::generalSettings().getExplorationOrder()), buildCommandLabels(false), buildAllRewardModels(false), buildStateInformation(false), rewardModelsToBuild(), constantDefinitions(), buildAllLabels(false), labelsToBuild(), expressionLabels(), terminalStates(), negatedTerminalStates() {
             if (formulas.empty()) {
                 this->buildAllRewardModels = true;
                 this->buildAllLabels = true;
@@ -253,6 +254,7 @@ namespace storm {
         template <typename ValueType, typename RewardModelType, typename StateType>
         std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::translate() {
             STORM_LOG_DEBUG("Building representation of program:" << std::endl << program << std::endl);
+            STORM_LOG_DEBUG("Exploration order is: " << options.explorationOrder);
             
             // Select the appropriate reward models (after the constants have been substituted).
             std::vector<std::reference_wrapper<storm::prism::RewardModel const>> selectedRewardModels;
@@ -314,7 +316,16 @@ namespace storm {
             std::pair<uint32_t, std::size_t> actualIndexBucketPair = internalStateInformation.stateStorage.findOrAddAndGetBucket(state, newIndex);
             
             if (actualIndexBucketPair.first == newIndex) {
-                statesToExplore.push_back(state);
+                if (options.explorationOrder == ExplorationOrder::Dfs) {
+                    statesToExplore.push_front(state);
+
+                    // Reserve one slot for the new state in the remapping.
+                    stateRemapping.get().push_back(storm::utility::zero<StateType>());
+                } else if (options.explorationOrder == ExplorationOrder::Bfs) {
+                    statesToExplore.push_back(state);
+                } else {
+                    STORM_LOG_ASSERT(false, "Invalid exploration order.");
+                }
                 ++internalStateInformation.numberOfStates;
             }
             
@@ -341,10 +352,18 @@ namespace storm {
             // Create a callback for the next-state generator to enable it to request the index of states.
             std::function<StateType (CompressedState const&)> stateToIdCallback = std::bind(&ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::getOrAddStateIndex, this, std::placeholders::_1);
             
+            // If the exploration order is something different from breadth-first, we need to keep track of the remapping
+            // from state ids to row groups. For this, we actually store the reversed mapping of row groups to state-ids
+            // and later reverse it.
+            if (options.explorationOrder != ExplorationOrder::Bfs) {
+                stateRemapping = std::vector<uint_fast64_t>();
+            }
+            
             // Let the generator create all initial states.
             this->internalStateInformation.initialStateIndices = generator.getInitialStates(stateToIdCallback);
             
             // Now explore the current state until there is no more reachable state.
+            uint_fast64_t currentRowGroup = 0;
             uint_fast64_t currentRow = 0;
 
             // Perform a search through the model.
@@ -354,6 +373,12 @@ namespace storm {
                 StateType currentIndex = internalStateInformation.stateStorage.getValue(currentState);
                 statesToExplore.pop_front();
                 
+                // If the exploration order differs from breadth-first, we remember that this row group was actually
+                // filled with the transitions of a different state.
+                if (options.explorationOrder != ExplorationOrder::Bfs) {
+                    stateRemapping.get()[currentIndex] = currentRowGroup;
+                }
+                
                 STORM_LOG_TRACE("Exploring state with id " << currentIndex << ".");
                 
                 storm::generator::StateBehavior<ValueType, StateType> behavior = generator.expand(currentState, stateToIdCallback);
@@ -384,6 +409,7 @@ namespace storm {
                         }
                         
                         ++currentRow;
+                        ++currentRowGroup;
                     } else {
                         std::cout << unpackStateIntoValuation(currentState).toString(true) << std::endl;
                         STORM_LOG_THROW(false, storm::exceptions::WrongFormatException,
@@ -430,9 +456,40 @@ namespace storm {
                         }
                         ++currentRow;
                     }
+                    ++currentRowGroup;
                 }
             }
 
+            // If the exploration order was not breadth-first, we need to fix the entries in the matrix according to
+            // (reversed) mapping of row groups to indices.
+            if (options.explorationOrder != ExplorationOrder::Bfs) {
+                STORM_LOG_ASSERT(stateRemapping, "Unable to fix columns without mapping.");
+                std::vector<uint_fast64_t> const& remapping = stateRemapping.get();
+                
+                // We need to fix the following entities:
+                // (a) the transition matrix
+                // (b) the initial states
+                // (c) the hash map storing the mapping states -> ids
+                
+                // Fix (a).
+                transitionMatrixBuilder.replaceColumns(remapping, 0);
+
+                // Fix (b).
+                std::vector<StateType> newInitialStateIndices(this->internalStateInformation.initialStateIndices.size());
+                std::transform(this->internalStateInformation.initialStateIndices.begin(), this->internalStateInformation.initialStateIndices.end(), newInitialStateIndices.begin(), [&remapping] (StateType const& state) { return remapping[state]; } );
+                std::sort(newInitialStateIndices.begin(), newInitialStateIndices.end());
+                this->internalStateInformation.initialStateIndices = std::move(newInitialStateIndices);
+                
+                // Fix (c).
+                std::unordered_map<StateType, StateType> valueRemapping;
+                for (StateType index = 0; index < remapping.size(); ++index) {
+                    if (remapping[index] != index) {
+                        valueRemapping.emplace(index, static_cast<StateType>(remapping[index]));
+                    }
+                }
+                this->internalStateInformation.stateStorage.remap(valueRemapping);
+            }
+            
             return choiceLabels;
         }
         
diff --git a/src/builder/ExplicitPrismModelBuilder.h b/src/builder/ExplicitPrismModelBuilder.h
index 3e5614854..ab2eb9815 100644
--- a/src/builder/ExplicitPrismModelBuilder.h
+++ b/src/builder/ExplicitPrismModelBuilder.h
@@ -1,5 +1,5 @@
-#ifndef STORM_ADAPTERS_EXPLICITPRISMMODELBUILDER_H
-#define	STORM_ADAPTERS_EXPLICITPRISMMODELBUILDER_H
+#ifndef STORM_BUILDER_EXPLICITPRISMMODELBUILDER_H
+#define	STORM_BUILDER_EXPLICITPRISMMODELBUILDER_H
 
 #include <memory>
 #include <utility>
@@ -24,6 +24,8 @@
 
 #include "src/utility/prism.h"
 
+#include "src/builder/ExplorationOrder.h"
+
 #include "src/generator/CompressedState.h"
 #include "src/generator/VariableInformation.h"
 
@@ -102,6 +104,11 @@ namespace storm {
                  */
                 Options();
                 
+                /*!
+                 * Copies the given set of options.
+                 */
+                Options(Options const& other) = default;
+                
                 /*! Creates an object representing the suggested building options assuming that the given formula is the
                  * only one to check. Additional formulas may be preserved by calling <code>preserveFormula</code>.
                  *
@@ -144,6 +151,9 @@ namespace storm {
                  */
                 void setTerminalStatesFromFormula(storm::logic::Formula const& formula);
                 
+                // The order in which to explore the model.
+                ExplorationOrder explorationOrder;
+                
                 // A flag that indicates whether or not command labels are to be built.
                 bool buildCommandLabels;
                 
@@ -286,13 +296,13 @@ namespace storm {
             // A set of states that still need to be explored.
             std::deque<CompressedState> statesToExplore;
             
-            // An optional mapping from row groups to the indices of the states that they reflect. This needs to be built
-            // in case the exploration order is not BFS.
-//            boost::optional<std::vector<StateType, StateType>> rowGroupToIndexMapping;
+            // An optional mapping from state indices to the row groups in which they actually reside. This needs to be
+            // built in case the exploration order is not BFS.
+            boost::optional<std::vector<uint_fast64_t>> stateRemapping;
 
         };
         
     } // namespace adapters
 } // namespace storm
 
-#endif	/* STORM_ADAPTERS_EXPLICITPRISMMODELBUILDER_H */
+#endif	/* STORM_BUILDER_EXPLICITPRISMMODELBUILDER_H */
diff --git a/src/builder/ExplorationOrder.cpp b/src/builder/ExplorationOrder.cpp
new file mode 100644
index 000000000..fb0b82766
--- /dev/null
+++ b/src/builder/ExplorationOrder.cpp
@@ -0,0 +1,22 @@
+#include "src/builder/ExplorationOrder.h"
+
+namespace storm {
+    namespace builder {
+        
+        std::ostream& operator<<(std::ostream& out, ExplorationOrder const& order) {
+            switch (order) {
+                case ExplorationOrder::Dfs:
+                    out << "depth-first";
+                    break;
+                case ExplorationOrder::Bfs:
+                    out << "breadth-first";
+                    break;
+                default:
+                    out << "undefined";
+                    break;
+            }
+            return out;
+        }
+        
+    }
+}
\ No newline at end of file
diff --git a/src/builder/ExplorationOrder.h b/src/builder/ExplorationOrder.h
new file mode 100644
index 000000000..94c6363b3
--- /dev/null
+++ b/src/builder/ExplorationOrder.h
@@ -0,0 +1,17 @@
+#ifndef STORM_BUILDER_EXPLORATIONORDER_H_
+#define STORM_BUILDER_EXPLORATIONORDER_H_
+
+#include <ostream>
+
+namespace storm {
+    namespace builder {
+        
+        // An enum that contains all currently supported exploration orders.
+        enum class ExplorationOrder { Dfs, Bfs };
+        
+        std::ostream& operator<<(std::ostream& out, ExplorationOrder const& order);
+        
+    }
+}
+
+#endif /* STORM_BUILDER_EXPLORATIONORDER_H_ */
\ No newline at end of file
diff --git a/src/settings/modules/GeneralSettings.cpp b/src/settings/modules/GeneralSettings.cpp
index 95a48931b..171b2a886 100644
--- a/src/settings/modules/GeneralSettings.cpp
+++ b/src/settings/modules/GeneralSettings.cpp
@@ -32,6 +32,8 @@ namespace storm {
             const std::string GeneralSettings::explicitOptionShortName = "exp";
             const std::string GeneralSettings::symbolicOptionName = "symbolic";
             const std::string GeneralSettings::symbolicOptionShortName = "s";
+            const std::string GeneralSettings::explorationOrderOptionName = "explorder";
+            const std::string GeneralSettings::explorationOrderOptionShortName = "eo";
             const std::string GeneralSettings::propertyOptionName = "prop";
             const std::string GeneralSettings::propertyOptionShortName = "prop";
             const std::string GeneralSettings::transitionRewardsOptionName = "transrew";
@@ -83,6 +85,11 @@ namespace storm {
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("labeling filename", "The name of the file from which to read the state labeling.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, symbolicOptionName, false, "Parses the model given in a symbolic representation.").setShortName(symbolicOptionShortName)
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The name of the file from which to read the symbolic model.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
+
+                std::vector<std::string> explorationOrders = {"dfs", "bfs"};
+                this->addOption(storm::settings::OptionBuilder(moduleName, explorationOrderOptionName, false, "Sets which exploration order to use.").setShortName(explorationOrderOptionShortName)
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the exploration order to choose. Available are: dfs and bfs.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(explorationOrders)).setDefaultValueString("dfs").build()).build());
+
                 this->addOption(storm::settings::OptionBuilder(moduleName, propertyOptionName, false, "Specifies the formulas to be checked on the model.").setShortName(propertyOptionShortName)
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("formula or filename", "The formula or the file containing the formulas.").build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, counterexampleOptionName, false, "Generates a counterexample for the given PRCTL formulas if not satisfied by the model")
@@ -186,6 +193,20 @@ namespace storm {
                 return this->getOption(symbolicOptionName).getArgumentByName("filename").getValueAsString();
             }
             
+            bool GeneralSettings::isExplorationOrderSet() const {
+                return this->getOption(explorationOrderOptionName).getHasOptionBeenSet();
+            }
+            
+            storm::builder::ExplorationOrder GeneralSettings::getExplorationOrder() const {
+                std::string explorationOrderAsString = this->getOption(explorationOrderOptionName).getArgumentByName("name").getValueAsString();
+                if (explorationOrderAsString == "dfs") {
+                    return storm::builder::ExplorationOrder::Dfs;
+                } else if (explorationOrderAsString == "bfs") {
+                    return storm::builder::ExplorationOrder::Bfs;
+                }
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown exploration order '" << explorationOrderAsString << "'.");
+            }
+            
             bool GeneralSettings::isPropertySet() const {
                 return this->getOption(propertyOptionName).getHasOptionBeenSet();
             }
diff --git a/src/settings/modules/GeneralSettings.h b/src/settings/modules/GeneralSettings.h
index 8144c0a1f..ad268db93 100644
--- a/src/settings/modules/GeneralSettings.h
+++ b/src/settings/modules/GeneralSettings.h
@@ -4,6 +4,8 @@
 #include "storm-config.h"
 #include "src/settings/modules/ModuleSettings.h"
 
+#include "src/builder/ExplorationOrder.h"
+
 namespace storm {
     namespace solver {
         enum class EquationSolverType;
@@ -25,7 +27,6 @@ namespace storm {
             class GeneralSettings : public ModuleSettings {
             public:
                 // An enumeration of all engines.
-
                 enum class Engine {
                     Sparse, Hybrid, Dd, AbstractionRefinement
                 };
@@ -138,6 +139,20 @@ namespace storm {
                  * @return The name of the file that contains the symbolic model specification.
                  */
                 std::string getSymbolicModelFilename() const;
+                
+                /*!
+                 * Retrieves whether the model exploration order was set.
+                 *
+                 * @return True if the model exploration option was set.
+                 */
+                bool isExplorationOrderSet() const;
+                
+                /*!
+                 * Retrieves the exploration order if it was set.
+                 *
+                 * @return The chosen exploration order.
+                 */
+                storm::builder::ExplorationOrder getExplorationOrder() const;
 
                 /*!
                  * Retrieves whether the property option was set.
@@ -389,6 +404,8 @@ namespace storm {
                 static const std::string explicitOptionShortName;
                 static const std::string symbolicOptionName;
                 static const std::string symbolicOptionShortName;
+                static const std::string explorationOrderOptionName;
+                static const std::string explorationOrderOptionShortName;
                 static const std::string propertyOptionName;
                 static const std::string propertyOptionShortName;
                 static const std::string transitionRewardsOptionName;
diff --git a/src/storage/BitVectorHashMap.cpp b/src/storage/BitVectorHashMap.cpp
index 881fc016a..fd4304eaa 100644
--- a/src/storage/BitVectorHashMap.cpp
+++ b/src/storage/BitVectorHashMap.cpp
@@ -233,6 +233,16 @@ namespace storm {
             return std::make_pair(buckets.get(bucket * bucketSize, bucketSize), values[bucket]);
         }
         
+        template<class ValueType, class Hash1, class Hash2>
+        void BitVectorHashMap<ValueType, Hash1, Hash2>::remap(std::unordered_map<ValueType, ValueType> const& remapping) {
+            for (auto pos : occupied) {
+                auto it = remapping.find(values[pos]);
+                if (it != remapping.end()) {
+                    values[pos] = it->second;
+                }
+            }
+        }
+        
         template class BitVectorHashMap<uint_fast64_t>;
         template class BitVectorHashMap<uint32_t>;
     }
diff --git a/src/storage/BitVectorHashMap.h b/src/storage/BitVectorHashMap.h
index c53928b09..7e8a1dc0e 100644
--- a/src/storage/BitVectorHashMap.h
+++ b/src/storage/BitVectorHashMap.h
@@ -3,6 +3,7 @@
 
 #include <cstdint>
 #include <functional>
+#include <unordered_map>
 
 #include "src/storage/BitVector.h"
 
@@ -123,6 +124,13 @@ namespace storm {
              */
             std::size_t capacity() const;
             
+            /*!
+             * Performs a remapping of all values stored by applying the given remapping.
+             *
+             * @param remapping The remapping to apply.
+             */
+            void remap(std::unordered_map<ValueType, ValueType> const& remapping);
+            
         private:
             /*!
              * Retrieves whether the given bucket holds a value.
diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp
index 3dbabf854..08138d48c 100644
--- a/src/storage/SparseMatrix.cpp
+++ b/src/storage/SparseMatrix.cpp
@@ -143,12 +143,12 @@ namespace storm {
             STORM_LOG_THROW(startingRow >= lastRow, storm::exceptions::InvalidStateException, "Illegal row group with negative size.");
             rowGroupIndices.push_back(startingRow);
             ++currentRowGroup;
-
+            
             // Close all rows from the most recent one to the starting row.
             for (index_type i = lastRow + 1; i <= startingRow; ++i) {
                 rowIndications.push_back(currentEntryCount);
             }
-
+            
             // Reset the most recently seen row/column to allow for proper insertion of the following elements.
             lastRow = startingRow;
             lastColumn = 0;
@@ -162,7 +162,7 @@ namespace storm {
                 rowCount = std::max(rowCount, initialRowCount);
             }
             rowCount = std::max(rowCount, overriddenRowCount);
-
+            
             // If the current row count was overridden, we may need to add empty rows.
             for (index_type i = lastRow + 1; i < rowCount; ++i) {
                 rowIndications.push_back(currentEntryCount);
@@ -202,7 +202,7 @@ namespace storm {
                     rowGroupIndices.push_back(rowCount);
                 }
             }
-
+            
             return SparseMatrix<ValueType>(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), hasCustomRowGrouping);
         }
         
@@ -216,6 +216,72 @@ namespace storm {
             return lastColumn;
         }
         
+        template<typename ValueType>
+        bool SparseMatrixBuilder<ValueType>::replaceColumns(std::vector<index_type> const& replacements, index_type offset) {
+            bool changed = false;
+            
+            // Sort columns per row.
+            typename SparseMatrix<ValueType>::index_type endGroups;
+            typename SparseMatrix<ValueType>::index_type endRows;
+            for (index_type group = 0; group < rowGroupIndices.size(); ++group) {
+                endGroups = group < rowGroupIndices.size()-1 ? rowGroupIndices[group+1] : rowIndications.size();
+                for (index_type i = rowGroupIndices[group]; i < endGroups; ++i) {
+                    bool changedRow = false;
+                    endRows = i < rowIndications.size()-1 ? rowIndications[i+1] : columnsAndValues.size();
+                    
+                    for (auto it = columnsAndValues.begin() + rowIndications[i], ite = columnsAndValues.begin() + endRows; it != ite; ++it) {
+                        if (it->getColumn() >= offset) {
+                            it->setColumn(replacements[it->getColumn() - offset]);
+                            changedRow = true;
+                        }
+                        // Update highest column in a way that only works if the highest appearing index does not become
+                        // lower during performing the replacement.
+                        highestColumn = std::max(highestColumn, it->getColumn());
+                    }
+                    
+                    if (changedRow) {
+                        changed = true;
+                        
+                        // Sort the row.
+                        std::sort(columnsAndValues.begin() + rowIndications[i], columnsAndValues.begin() + endRows,
+                                  [](MatrixEntry<index_type, value_type> const& a, MatrixEntry<index_type, value_type> const& b) {
+                                      return a.getColumn() < b.getColumn();
+                                  });
+                        // Assert no equal elements
+                        STORM_LOG_ASSERT(std::is_sorted(columnsAndValues.begin() + rowIndications[i], columnsAndValues.begin() + endRows,
+                                                        [](MatrixEntry<index_type, value_type> const& a, MatrixEntry<index_type, value_type> const& b) {
+                                                            return a.getColumn() <= b.getColumn();
+                                                        }), "Must not have different elements with the same column in a row.");
+                    }
+                }
+            }
+            
+            return changed;
+        }
+        
+        template<typename ValueType>
+        void SparseMatrixBuilder<ValueType>::fixColumns() {
+            // Sort columns per row.
+            typename SparseMatrix<ValueType>::index_type endGroups;
+            typename SparseMatrix<ValueType>::index_type endRows;
+            for (index_type group = 0; group < rowGroupIndices.size(); ++group) {
+                endGroups = group < rowGroupIndices.size()-1 ? rowGroupIndices[group+1] : rowIndications.size();
+                for (index_type i = rowGroupIndices[group]; i < endGroups; ++i) {
+                    endRows = i < rowIndications.size()-1 ? rowIndications[i+1] : columnsAndValues.size();
+                    // Sort the row.
+                    std::sort(columnsAndValues.begin() + rowIndications[i], columnsAndValues.begin() + endRows,
+                              [](MatrixEntry<index_type, value_type> const& a, MatrixEntry<index_type, value_type> const& b) {
+                                  return a.getColumn() < b.getColumn();
+                              });
+                    // Assert no equal elements
+                    STORM_LOG_ASSERT(std::is_sorted(columnsAndValues.begin() + rowIndications[i], columnsAndValues.begin() + endRows,
+                                                    [](MatrixEntry<index_type, value_type> const& a, MatrixEntry<index_type, value_type> const& b) {
+                                                        return a.getColumn() <= b.getColumn();
+                                                    }), "Must not have different elements with the same column in a row.");
+                }
+            }
+        }
+        
         template<typename ValueType>
         SparseMatrix<ValueType>::rows::rows(iterator begin, index_type entryCount) : beginIterator(begin), entryCount(entryCount) {
             // Intentionally left empty.
@@ -532,12 +598,12 @@ namespace storm {
                     *it = i;
                 }
                 auto res = getSubmatrix(rowConstraint, columnConstraint, fakeRowGroupIndices, insertDiagonalElements);
-
+                
                 // Create a new row grouping that reflects the new sizes of the row groups.
                 std::vector<uint_fast64_t> newRowGroupIndices;
                 newRowGroupIndices.push_back(0);
                 auto selectedRowIt = rowConstraint.begin();
-
+                
                 // For this, we need to count how many rows were preserved in every group.
                 for (uint_fast64_t group = 0; group < this->getRowGroupCount(); ++group) {
                     uint_fast64_t newRowCount = 0;
@@ -1172,7 +1238,7 @@ namespace storm {
                 assert(this->getRowGroupSize(group) == 1);
                 for (typename SparseMatrix<ValueType>::index_type i = this->getRowGroupIndices()[group]; i < this->getRowGroupIndices()[group + 1]; ++i) {
                     typename SparseMatrix<ValueType>::index_type nextIndex = this->rowIndications[i];
-
+                    
                     // Print the actual row.
                     out << i << "\t(";
                     typename SparseMatrix<ValueType>::index_type currentRealIndex = 0;
@@ -1226,7 +1292,7 @@ namespace storm {
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<double> const& matrix);
         template std::vector<double> SparseMatrix<double>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<double> const& otherMatrix) const;
         template bool SparseMatrix<double>::isSubmatrixOf(SparseMatrix<double> const& matrix) const;
-
+        
         // float
         template class MatrixEntry<typename SparseMatrix<float>::index_type, float>;
         template std::ostream& operator<<(std::ostream& out, MatrixEntry<typename SparseMatrix<float>::index_type, float> const& entry);
@@ -1235,7 +1301,7 @@ namespace storm {
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<float> const& matrix);
         template std::vector<float> SparseMatrix<float>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<float> const& otherMatrix) const;
         template bool SparseMatrix<float>::isSubmatrixOf(SparseMatrix<float> const& matrix) const;
-
+        
         // int
         template class MatrixEntry<typename SparseMatrix<int>::index_type, int>;
         template std::ostream& operator<<(std::ostream& out, MatrixEntry<typename SparseMatrix<int>::index_type, int> const& entry);
@@ -1243,7 +1309,7 @@ namespace storm {
         template class SparseMatrix<int>;
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<int> const& matrix);
         template bool SparseMatrix<int>::isSubmatrixOf(SparseMatrix<int> const& matrix) const;
-
+        
         // state_type
         template class MatrixEntry<typename SparseMatrix<storm::storage::sparse::state_type>::index_type, storm::storage::sparse::state_type>;
         template std::ostream& operator<<(std::ostream& out, MatrixEntry<typename SparseMatrix<storm::storage::sparse::state_type>::index_type, storm::storage::sparse::state_type> const& entry);
@@ -1251,7 +1317,7 @@ namespace storm {
         template class SparseMatrix<storm::storage::sparse::state_type>;
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<storm::storage::sparse::state_type> const& matrix);
         template bool SparseMatrix<int>::isSubmatrixOf(SparseMatrix<storm::storage::sparse::state_type> const& matrix) const;
-
+        
 #ifdef STORM_HAVE_CARL
         // Rat Function
         template class MatrixEntry<typename SparseMatrix<RationalFunction>::index_type, RationalFunction>;
@@ -1264,7 +1330,7 @@ namespace storm {
         template std::vector<storm::RationalFunction> SparseMatrix<float>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<storm::RationalFunction> const& otherMatrix) const;
         template std::vector<storm::RationalFunction> SparseMatrix<int>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<storm::RationalFunction> const& otherMatrix) const;
         template bool SparseMatrix<storm::RationalFunction>::isSubmatrixOf(SparseMatrix<storm::RationalFunction> const& matrix) const;
-
+        
         // Intervals
         template std::vector<storm::Interval> SparseMatrix<double>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<storm::Interval> const& otherMatrix) const;
         template class MatrixEntry<typename SparseMatrix<Interval>::index_type, Interval>;
@@ -1274,7 +1340,7 @@ namespace storm {
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<Interval> const& matrix);
         template std::vector<storm::Interval> SparseMatrix<Interval>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<storm::Interval> const& otherMatrix) const;
         template bool SparseMatrix<storm::Interval>::isSubmatrixOf(SparseMatrix<storm::Interval> const& matrix) const;
-
+        
         template bool SparseMatrix<storm::Interval>::isSubmatrixOf(SparseMatrix<double> const& matrix) const;
 #endif
         
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index f9d9916c4..4892dbe32 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -217,6 +217,17 @@ namespace storm {
              */
             index_type getLastColumn() const;
             
+            /*!
+             * Replaces all columns with id > offset according to replacements.
+             * Every state  with id offset+i is replaced by the id in replacements[i].
+             * Afterwards the columns are sorted.
+             *
+             * @param replacements Mapping indicating the replacements from offset+i -> value of i.
+             * @param offset Offset to add to each id in vector index.
+             * @return True if replacement took place, False if nothing changed.
+             */
+            bool replaceColumns(std::vector<index_type> const& replacements, index_type offset);
+            
         private:
             // A flag indicating whether a row count was set upon construction.
             bool initialRowCountSet;
@@ -277,6 +288,11 @@ namespace storm {
             // Stores the currently active row group. This is used for correctly constructing the row grouping of the
             // matrix.
             index_type currentRowGroup;
+            
+            /*!
+             * Fixes the matrix by sorting the columns to gain increasing order again.
+             */
+            void fixColumns();
         };
         
         /*!
diff --git a/src/storage/prism/RewardModel.cpp b/src/storage/prism/RewardModel.cpp
index 3f38f7bd8..30c43728c 100644
--- a/src/storage/prism/RewardModel.cpp
+++ b/src/storage/prism/RewardModel.cpp
@@ -84,7 +84,7 @@ namespace storm {
         std::ostream& operator<<(std::ostream& stream, RewardModel const& rewardModel) {
             stream << "rewards";
             if (rewardModel.getName() != "") {
-                std::cout << " \"" << rewardModel.getName() << "\"";
+                stream << " \"" << rewardModel.getName() << "\"";
             }
             stream << std::endl;
             for (auto const& reward : rewardModel.getStateRewards()) {
diff --git a/src/utility/constants.cpp b/src/utility/constants.cpp
index 900482a72..941f07f7e 100644
--- a/src/utility/constants.cpp
+++ b/src/utility/constants.cpp
@@ -192,6 +192,10 @@ namespace storm {
         template bool isZero(storm::storage::sparse::state_type const& value);
         template bool isConstant(storm::storage::sparse::state_type const& value);
 
+        template uint32_t one();
+        template uint32_t zero();
+        template uint32_t infinity();
+                
         template storm::storage::sparse::state_type one();
         template storm::storage::sparse::state_type zero();
         template storm::storage::sparse::state_type infinity();