allow arbitrary-layer meta variables in DdManager; make DdManager available as non-const from a DD; started on symbolic state elimination linear equation solver

src/storm/models/symbolic/Model.cpp

@@ -13,6 +13,8 @@
 
 #include "storm/models/symbolic/StandardRewardModel.h"
 
+#include "storm/utility/dd.h"
+
 #include "storm-config.h"
 #include "storm/adapters/CarlAdapter.h"
 
@@ -146,12 +148,7 @@ namespace storm {
             
             template<storm::dd::DdType Type, typename ValueType>
             storm::dd::Add<Type, ValueType> Model<Type, ValueType>::getRowColumnIdentity() const {
-                storm::dd::Add<Type, ValueType> result = this->getManager().template getAddOne<ValueType>();
-                for (auto const& pair : this->getRowColumnMetaVariablePairs()) {
-                    result *= this->getManager().template getIdentity<ValueType>(pair.first).equals(this->getManager().template getIdentity<ValueType>(pair.second)).template toAdd<ValueType>();
-                    result *= this->getManager().getRange(pair.first).template toAdd<ValueType>() * this->getManager().getRange(pair.second).template toAdd<ValueType>();
-                }
-                return result;
+                return storm::utility::dd::getRowColumnDiagonal<Type, ValueType>(this->getManager(), this->getRowColumnMetaVariablePairs());
             }
             
             template<storm::dd::DdType Type, typename ValueType>

src/storm/solver/EliminationLinearEquationSolver.cpp

@@ -109,7 +109,7 @@ namespace storm {
             // After having solved the system, we need to revert the transition system if we kept it local.
             if (localA) {
                 localA->convertToEquationSystem();
-            };
+            }
 
             return true;
         }

src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp

@@ -0,0 +1,64 @@
+#include "storm/solver/SymbolicEliminationLinearEquationSolver.h"
+
+#include "storm/storage/dd/DdManager.h"
+#include "storm/storage/dd/Add.h"
+
+#include "storm/utility/dd.h"
+
+namespace storm {
+    namespace solver {
+     
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicEliminationLinearEquationSolver<DdType, ValueType>::SymbolicEliminationLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) : SymbolicLinearEquationSolver<DdType, ValueType>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicEliminationLinearEquationSolver<DdType, ValueType>::SymbolicEliminationLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : SymbolicLinearEquationSolver<DdType, ValueType>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, precision, maximalNumberOfIterations, relative) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType, ValueType> SymbolicEliminationLinearEquationSolver<DdType, ValueType>::solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
+            storm::dd::Bdd<DdType> diagonal = storm::utility::dd::getRowColumnDiagonal(x.getDdManager(), this->rowColumnMetaVariablePairs);
+            diagonal &= this->allRows;
+            
+            storm::dd::Add<DdType, ValueType> rowsAdd = this->allRows.template toAdd<ValueType>();
+            storm::dd::Add<DdType, ValueType> diagonalAdd = diagonal.template toAdd<ValueType>();
+            
+            // Revert the conversion to an equation system.
+            storm::dd::Add<DdType, ValueType> matrix = diagonalAdd - this->A;
+            
+            storm::dd::Add<DdType, ValueType> solution = b;
+            
+            // As long as there are transitions, we eliminate them.
+            while (!matrix.isZero()) {
+                // Determine inverse loop probabilies
+                storm::dd::Add<DdType, ValueType> inverseLoopProbabilities = rowsAdd / (rowsAdd - (diagonalAdd * matrix).sumAbstract(this->columnMetaVariables));
+                
+                inverseLoopProbabilities.swapVariables(this->rowColumnMetaVariablePairs);
+                
+                // Scale all transitions with the inverse loop probabilities.
+                matrix *= inverseLoopProbabilities;
+                
+                // Delete diagonal elements, i.e. remove self-loops.
+                matrix = diagonal.ite(x.getDdManager().template getAddZero<ValueType>(), matrix);
+            
+                // Update the one-step probabilities.
+                solution += (matrix * solution.swapVariables(this->rowColumnMetaVariablePairs)).sumAbstract(this->columnMetaVariables);
+                
+                // Now eliminate all direct transitions of all states.
+                storm::dd::Add<DdType, ValueType> matrixWithRemoved;
+            }
+            
+            std::cout << "here" << std::endl;
+            solution.exportToDot("solution.dot");
+            
+            exit(-1);
+        }
+
+        template class SymbolicEliminationLinearEquationSolver<storm::dd::DdType::CUDD, double>;
+        template class SymbolicEliminationLinearEquationSolver<storm::dd::DdType::Sylvan, double>;
+        
+    }
+}

src/storm/solver/SymbolicLinearEquationSolver.cpp

@@ -3,6 +3,8 @@
 #include "storm/storage/dd/DdManager.h"
 #include "storm/storage/dd/Add.h"
 
+#include "storm/utility/dd.h"
+
 #include "storm/settings/SettingsManager.h"
 #include "storm/settings/modules/NativeEquationSolverSettings.h"
 
@@ -28,11 +30,7 @@ namespace storm {
         template<storm::dd::DdType DdType, typename ValueType>
         storm::dd::Add<DdType, ValueType>  SymbolicLinearEquationSolver<DdType, ValueType>::solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
             // Start by computing the Jacobi decomposition of the matrix A.
-            storm::dd::Bdd<DdType> diagonal = x.getDdManager().getBddOne();
-            for (auto const& pair : rowColumnMetaVariablePairs) {
-                diagonal &= x.getDdManager().template getIdentity<ValueType>(pair.first).equals(x.getDdManager().template getIdentity<ValueType>(pair.second));
-                diagonal &= x.getDdManager().getRange(pair.first) && x.getDdManager().getRange(pair.second);
-            }
+            storm::dd::Bdd<DdType> diagonal = storm::utility::dd::getRowColumnDiagonal(x.getDdManager(), rowColumnMetaVariablePairs);
             diagonal &= allRows;
             
             storm::dd::Add<DdType, ValueType> lu = diagonal.ite(this->A.getDdManager().template getAddZero<ValueType>(), this->A);

src/storm/solver/SymbolicLinearEquationSolver.h

@@ -1,15 +1,12 @@
 #ifndef STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_
 #define STORM_SOLVER_SYMBOLICLINEAREQUATIONSOLVER_H_
 
-#include <memory>
 #include <set>
 #include <vector>
-#include <boost/variant.hpp>
 
 #include "storm/storage/expressions/Variable.h"
 #include "storm/storage/dd/DdType.h"
 
-
 namespace storm {
     namespace dd {
         template<storm::dd::DdType Type, typename ValueType>

src/storm/storage/dd/Dd.cpp

@@ -10,7 +10,7 @@
 namespace storm {
     namespace dd {
         template<DdType LibraryType>
-        Dd<LibraryType>::Dd(DdManager<LibraryType> const& ddManager, std::set<storm::expressions::Variable> const& containedMetaVariables) : ddManager(&ddManager), containedMetaVariables(containedMetaVariables) {
+        Dd<LibraryType>::Dd(DdManager<LibraryType> const& ddManager, std::set<storm::expressions::Variable> const& containedMetaVariables) : ddManager(const_cast<DdManager<LibraryType>*>(&ddManager)), containedMetaVariables(containedMetaVariables) {
             // Intentionally left empty.
         }
         
@@ -35,7 +35,7 @@ namespace storm {
         }
         
         template<DdType LibraryType>
-        DdManager<LibraryType> const& Dd<LibraryType>::getDdManager() const {
+        DdManager<LibraryType>& Dd<LibraryType>::getDdManager() const {
             return *this->ddManager;
         }
         

src/storm/storage/dd/Dd.h

@@ -109,7 +109,7 @@ namespace storm {
              *
              * A pointer to the manager that is responsible for this DD.
              */
-            DdManager<LibraryType> const& getDdManager() const;
+            DdManager<LibraryType>& getDdManager() const;
             
             /*!
              * Retrieves the set of all meta variables contained in the DD.
@@ -182,7 +182,7 @@ namespace storm {
             
         private:
             // A pointer to the manager responsible for this DD.
-            DdManager<LibraryType> const* ddManager;
+            DdManager<LibraryType>* ddManager;
             
             // The meta variables that appear in this DD.
             std::set<storm::expressions::Variable> containedMetaVariables;

src/storm/storage/dd/DdManager.cpp

@@ -126,6 +126,15 @@ namespace storm {
         
         template<DdType LibraryType>
         std::pair<storm::expressions::Variable, storm::expressions::Variable> DdManager<LibraryType>::addMetaVariable(std::string const& name, int_fast64_t low, int_fast64_t high, boost::optional<std::pair<MetaVariablePosition, storm::expressions::Variable>> const& position) {
+            std::vector<storm::expressions::Variable> result = addMetaVariable(name, low, high, 2, position);
+            return std::make_pair(result[0], result[1]);
+        }
+        
+        template<DdType LibraryType>
+        std::vector<storm::expressions::Variable> DdManager<LibraryType>::addMetaVariable(std::string const& name, int_fast64_t low, int_fast64_t high, uint64_t numberOfLayers, boost::optional<std::pair<MetaVariablePosition, storm::expressions::Variable>> const& position) {
+            // Check whether number of layers is legal.
+            STORM_LOG_THROW(numberOfLayers >= 1, storm::exceptions::InvalidArgumentException, "Layers must be at least 1.");
+            
             // Check whether the variable name is legal.
             STORM_LOG_THROW(name != "" && name.back() != '\'', storm::exceptions::InvalidArgumentException, "Illegal name of meta variable: '" << name << "'.");
             
@@ -155,31 +164,48 @@ namespace storm {
                 ++numberOfBits;
             }
             
-            storm::expressions::Variable unprimed = manager->declareBitVectorVariable(name, numberOfBits);
-            storm::expressions::Variable primed = manager->declareBitVectorVariable(name + "'", numberOfBits);
+            std::stringstream tmp1;
+            std::vector<storm::expressions::Variable> result;
+            for (uint64 layer = 0; layer < numberOfLayers; ++layer) {
+                result.emplace_back(manager->declareBitVectorVariable(name + tmp1.str(), numberOfBits));
+                tmp1 << "'";
+            }
+            
+            std::vector<std::vector<Bdd<LibraryType>>> variables(numberOfLayers);
             
-            std::vector<Bdd<LibraryType>> variables;
-            std::vector<Bdd<LibraryType>> variablesPrime;
             for (std::size_t i = 0; i < numberOfBits; ++i) {
-                auto ddVariablePair = internalDdManager.createNewDdVariablePair(level);
-                variables.emplace_back(Bdd<LibraryType>(*this, ddVariablePair.first, {unprimed}));
-                variablesPrime.emplace_back(Bdd<LibraryType>(*this, ddVariablePair.second, {primed}));
+                std::vector<InternalBdd<LibraryType>> ddVariables = internalDdManager.createDdVariables(numberOfLayers, level);
+                for (uint64 layer = 0; layer < numberOfLayers; ++layer) {
+                    variables[layer].emplace_back(Bdd<LibraryType>(*this, ddVariables[layer], {result[layer]}));
+                }
                 
                 // If we are inserting the variable at a specific level, we need to prepare the level for the next pair
                 // of variables.
                 if (level) {
-                    level.get() += 2;
+                    level.get() += numberOfLayers;
                 }
             }
-
-            metaVariableMap.emplace(unprimed, DdMetaVariable<LibraryType>(name, low, high, variables));
-            metaVariableMap.emplace(primed, DdMetaVariable<LibraryType>(name + "'", low, high, variablesPrime));
             
-            return std::make_pair(unprimed, primed);
+            std::stringstream tmp2;
+            for (uint64_t layer = 0; layer < numberOfLayers; ++layer) {
+                metaVariableMap.emplace(result[layer], DdMetaVariable<LibraryType>(name + tmp2.str(), low, high, variables[layer]));
+                tmp2 << "'";
+            }
+            
+            return result;
         }
         
         template<DdType LibraryType>
         std::pair<storm::expressions::Variable, storm::expressions::Variable> DdManager<LibraryType>::addMetaVariable(std::string const& name, boost::optional<std::pair<MetaVariablePosition, storm::expressions::Variable>> const& position) {
+            std::vector<storm::expressions::Variable> result = addMetaVariable(name, 2, position);
+            return std::make_pair(result[0], result[1]);
+        }
+        
+        template<DdType LibraryType>
+        std::vector<storm::expressions::Variable> DdManager<LibraryType>::addMetaVariable(std::string const& name, uint64_t numberOfLayers, boost::optional<std::pair<MetaVariablePosition, storm::expressions::Variable>> const& position) {
+            // Check whether number of layers is legal.
+            STORM_LOG_THROW(numberOfLayers >= 1, storm::exceptions::InvalidArgumentException, "Layers must be at least 1.");
+
             // Check whether the variable name is legal.
             STORM_LOG_THROW(name != "" && name.back() != '\'', storm::exceptions::InvalidArgumentException, "Illegal name of meta variable: '" << name << "'.");
             
@@ -200,19 +226,28 @@ namespace storm {
                 }
             }
             
-            storm::expressions::Variable unprimed = manager->declareBooleanVariable(name);
-            storm::expressions::Variable primed = manager->declareBooleanVariable(name + "'");
+            std::stringstream tmp1;
+            std::vector<storm::expressions::Variable> result;
+            for (uint64 layer = 0; layer < numberOfLayers; ++layer) {
+                result.emplace_back(manager->declareBooleanVariable(name + tmp1.str()));
+                tmp1 << "'";
+            }
             
-            std::vector<Bdd<LibraryType>> variables;
-            std::vector<Bdd<LibraryType>> variablesPrime;
-            auto ddVariablePair = internalDdManager.createNewDdVariablePair(level);
-            variables.emplace_back(Bdd<LibraryType>(*this, ddVariablePair.first, {unprimed}));
-            variablesPrime.emplace_back(Bdd<LibraryType>(*this, ddVariablePair.second, {primed}));
+            std::vector<std::vector<Bdd<LibraryType>>> variables(numberOfLayers);
             
-            metaVariableMap.emplace(unprimed, DdMetaVariable<LibraryType>(name, variables));
-            metaVariableMap.emplace(primed, DdMetaVariable<LibraryType>(name + "'", variablesPrime));
+            std::vector<InternalBdd<LibraryType>> ddVariables = internalDdManager.createDdVariables(numberOfLayers, level);
             
-            return std::make_pair(unprimed, primed);
+            for (uint64_t layer = 0; layer < numberOfLayers; ++layer) {
+                variables[layer].emplace_back(Bdd<LibraryType>(*this, ddVariables[layer], {result[layer]}));
+            }
+            
+            std::stringstream tmp2;
+            for (uint64_t layer = 0; layer < numberOfLayers; ++layer) {
+                metaVariableMap.emplace(result[layer], DdMetaVariable<LibraryType>(name + tmp2.str(), variables[layer]));
+                tmp2 << "'";
+            }
+
+            return result;
         }
         
         template<DdType LibraryType>

src/storm/storage/dd/DdManager.h

@@ -135,21 +135,37 @@ namespace storm {
             Bdd<LibraryType> getCube(std::set<storm::expressions::Variable> const& variables) const;
 
             /*!
-             * Adds an integer meta variable with the given range.
+             * Adds an integer meta variable with the given range with two layers (a 'normal' and a 'primed' one).
              *
              * @param variableName The name of the new variable.
              * @param low The lowest value of the range of the variable.
              * @param high The highest value of the range of the variable.
              */
             std::pair<storm::expressions::Variable, storm::expressions::Variable> addMetaVariable(std::string const& variableName, int_fast64_t low, int_fast64_t high, boost::optional<std::pair<MetaVariablePosition, storm::expressions::Variable>> const& position = boost::none);
+
+            /*!
+             * Creates a meta variable with the given number of layers.
+             *
+             * @param variableName The name of the variable.
+             * @param numberOfLayers The number of layers of this variable (must be greater or equal 1).
+             */
+            std::vector<storm::expressions::Variable> addMetaVariable(std::string const& variableName, int_fast64_t low, int_fast64_t high, uint64_t numberOfLayers, boost::optional<std::pair<MetaVariablePosition, storm::expressions::Variable>> const& position = boost::none);
             
             /*!
-             * Adds a boolean meta variable.
+             * Adds a boolean meta variable with two layers (a 'normal' and a 'primed' one).
              *
              * @param variableName The name of the new variable.
              */
             std::pair<storm::expressions::Variable, storm::expressions::Variable> addMetaVariable(std::string const& variableName, boost::optional<std::pair<MetaVariablePosition, storm::expressions::Variable>> const& position = boost::none);
             
+            /*!
+             * Creates a meta variable with the given number of layers.
+             *
+             * @param variableName The name of the variable.
+             * @param numberOfLayers The number of layers of this variable (must be greater or equal 1).
+             */
+            std::vector<storm::expressions::Variable> addMetaVariable(std::string const& variableName, uint64_t numberOfLayers, boost::optional<std::pair<MetaVariablePosition, storm::expressions::Variable>> const& position = boost::none);
+            
             /*!
              * Retrieves the names of all meta variables that have been added to the manager.
              *

src/storm/storage/dd/cudd/InternalCuddDdManager.cpp

@@ -61,22 +61,24 @@ namespace storm {
             return InternalAdd<DdType::CUDD, ValueType>(this, cuddManager.constant(value));
         }
 
-        std::pair<InternalBdd<DdType::CUDD>, InternalBdd<DdType::CUDD>> InternalDdManager<DdType::CUDD>::createNewDdVariablePair(boost::optional<uint_fast64_t> const& position) {
-            std::pair<InternalBdd<DdType::CUDD>, InternalBdd<DdType::CUDD>> result;
+        std::vector<InternalBdd<DdType::CUDD>> InternalDdManager<DdType::CUDD>::createDdVariables(uint64_t numberOfLayers, boost::optional<uint_fast64_t> const& position) {
+            std::vector<InternalBdd<DdType::CUDD>> result;
             
             if (position) {
-                result.first = InternalBdd<DdType::CUDD>(this, cuddManager.bddNewVarAtLevel(position.get()));
-                result.second = InternalBdd<DdType::CUDD>(this, cuddManager.bddNewVarAtLevel(position.get() + 1));
+                for (uint64_t layer = 0; layer < numberOfLayers; ++layer) {
+                    result.emplace_back(InternalBdd<DdType::CUDD>(this, cuddManager.bddNewVarAtLevel(position.get() + layer)));
+                }
             } else {
-                result.first = InternalBdd<DdType::CUDD>(this, cuddManager.bddVar());
-                result.second = InternalBdd<DdType::CUDD>(this, cuddManager.bddVar());
+                for (uint64_t layer = 0; layer < numberOfLayers; ++layer) {
+                    result.emplace_back(InternalBdd<DdType::CUDD>(this, cuddManager.bddVar()));
+                }
             }
             
             // Connect the two variables so they are not 'torn apart' during dynamic reordering.
-            cuddManager.MakeTreeNode(result.first.getIndex(), 2, MTR_FIXED);
+            cuddManager.MakeTreeNode(result.front().getIndex(), numberOfLayers, MTR_FIXED);
             
             // Keep track of the number of variables.
-            numberOfDdVariables += 2;
+            numberOfDdVariables += numberOfLayers;
             
             return result;
         }

src/storm/storage/dd/cudd/InternalCuddDdManager.h

@@ -76,13 +76,13 @@ namespace storm {
             InternalAdd<DdType::CUDD, ValueType> getConstant(ValueType const& value) const;
             
             /*!
-             * Creates a new pair of DD variables and returns the two cubes as a result.
+             * Creates new layered DD variables and returns the cubes as a result.
              *
-             * @param position An optional position at which to insert the new variable pair. This may only be given, if
-             * the manager supports ordered insertion.
-             * @return The two cubes belonging to the DD variables.
+             * @param position An optional position at which to insert the new variable. This may only be given, if the
+             * manager supports ordered insertion.
+             * @return The cubes belonging to the DD variables.
              */
-            std::pair<InternalBdd<DdType::CUDD>, InternalBdd<DdType::CUDD>> createNewDdVariablePair(boost::optional<uint_fast64_t> const& position = boost::none);
+            std::vector<InternalBdd<DdType::CUDD>> createDdVariables(uint64_t numberOfLayers, boost::optional<uint_fast64_t> const& position = boost::none);
             
             /*!
              * Checks whether this manager supports the ordered insertion of variables, i.e. inserting variables at

src/storm/storage/dd/sylvan/InternalSylvanDdManager.cpp

@@ -126,13 +126,17 @@ namespace storm {
 		}
 #endif
 
-        std::pair<InternalBdd<DdType::Sylvan>, InternalBdd<DdType::Sylvan>> InternalDdManager<DdType::Sylvan>::createNewDdVariablePair(boost::optional<uint_fast64_t> const& position) {
+        std::vector<InternalBdd<DdType::Sylvan>> InternalDdManager<DdType::Sylvan>::createDdVariables(uint64_t numberOfLayers, boost::optional<uint_fast64_t> const& position) {
             STORM_LOG_THROW(!position, storm::exceptions::NotSupportedException, "The manager does not support ordered insertion.");
 			
-            InternalBdd<DdType::Sylvan> first = InternalBdd<DdType::Sylvan>(this, sylvan::Bdd::bddVar(nextFreeVariableIndex));
-            InternalBdd<DdType::Sylvan> second = InternalBdd<DdType::Sylvan>(this, sylvan::Bdd::bddVar(nextFreeVariableIndex + 1));
-            nextFreeVariableIndex += 2;
-            return std::make_pair(first, second);
+            std::vector<InternalBdd<DdType::Sylvan>> result;
+            
+            for (uint64_t layer = 0; layer < numberOfLayers; ++layer) {
+                result.emplace_back(InternalBdd<DdType::Sylvan>(this, sylvan::Bdd::bddVar(nextFreeVariableIndex)));
+                ++nextFreeVariableIndex;
+            }
+            
+            return result;
         }
         
         bool InternalDdManager<DdType::Sylvan>::supportsOrderedInsertion() const {

src/storm/storage/dd/sylvan/InternalSylvanDdManager.h

@@ -77,13 +77,13 @@ namespace storm {
             InternalAdd<DdType::Sylvan, ValueType> getConstant(ValueType const& value) const;
             
             /*!
-             * Creates a new pair of DD variables and returns the two cubes as a result.
+             * Creates new layered DD variables and returns the cubes as a result.
              *
-             * @param position An optional position at which to insert the new variable pair. This may only be given, if
-             * the manager supports ordered insertion.
-             * @return The two cubes belonging to the DD variables.
+             * @param position An optional position at which to insert the new variable. This may only be given, if the
+             * manager supports ordered insertion.
+             * @return The cubes belonging to the DD variables.
              */
-            std::pair<InternalBdd<DdType::Sylvan>, InternalBdd<DdType::Sylvan>> createNewDdVariablePair(boost::optional<uint_fast64_t> const& position = boost::none);
+            std::vector<InternalBdd<DdType::Sylvan>> createDdVariables(uint64_t numberOfLayers, boost::optional<uint_fast64_t> const& position = boost::none);
             
             /*!
              * Checks whether this manager supports the ordered insertion of variables, i.e. inserting variables at

src/storm/utility/dd.cpp

@@ -1,8 +1,11 @@
 #include "storm/utility/dd.h"
 
+#include "storm/storage/dd/DdManager.h"
 #include "storm/storage/dd/Add.h"
 #include "storm/storage/dd/Bdd.h"
 
+#include "storm/adapters/CarlAdapter.h"
+
 #include "storm/utility/macros.h"
 
 namespace storm {
@@ -42,9 +45,37 @@ namespace storm {
                 return reachableStates;
             }
             
+            template <storm::dd::DdType Type, typename ValueType>
+            storm::dd::Add<Type, ValueType> getRowColumnDiagonal(storm::dd::DdManager<Type> const& ddManager, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) {
+                storm::dd::Add<Type, ValueType> result = ddManager.template getAddOne<ValueType>();
+                for (auto const& pair : rowColumnMetaVariablePairs) {
+                    result *= ddManager.template getIdentity<ValueType>(pair.first).equals(ddManager.template getIdentity<ValueType>(pair.second)).template toAdd<ValueType>();
+                    result *= ddManager.getRange(pair.first).template toAdd<ValueType>() * ddManager.getRange(pair.second).template toAdd<ValueType>();
+                }
+                return result;
+            }
+            
+            template <storm::dd::DdType Type>
+            storm::dd::Bdd<Type> getRowColumnDiagonal(storm::dd::DdManager<Type> const& ddManager, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) {
+                storm::dd::Bdd<Type> diagonal = ddManager.getBddOne();
+                for (auto const& pair : rowColumnMetaVariablePairs) {
+                    diagonal &= ddManager.template getIdentity<uint64_t>(pair.first).equals(ddManager.template getIdentity<uint64_t>(pair.second));
+                    diagonal &= ddManager.getRange(pair.first) && ddManager.getRange(pair.second);
+                }
+                return diagonal;
+            }
+            
             template storm::dd::Bdd<storm::dd::DdType::CUDD> computeReachableStates(storm::dd::Bdd<storm::dd::DdType::CUDD> const& initialStates, storm::dd::Bdd<storm::dd::DdType::CUDD> const& transitions, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables);
             template storm::dd::Bdd<storm::dd::DdType::Sylvan> computeReachableStates(storm::dd::Bdd<storm::dd::DdType::Sylvan> const& initialStates, storm::dd::Bdd<storm::dd::DdType::Sylvan> const& transitions, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables);
-            
+
+            template storm::dd::Add<storm::dd::DdType::CUDD, double> getRowColumnDiagonal(storm::dd::DdManager<storm::dd::DdType::CUDD> const& ddManager, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+            template storm::dd::Add<storm::dd::DdType::Sylvan, double> getRowColumnDiagonal(storm::dd::DdManager<storm::dd::DdType::Sylvan> const& ddManager, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+
+            template storm::dd::Bdd<storm::dd::DdType::CUDD> getRowColumnDiagonal(storm::dd::DdManager<storm::dd::DdType::CUDD> const& ddManager, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+            template storm::dd::Bdd<storm::dd::DdType::Sylvan> getRowColumnDiagonal(storm::dd::DdManager<storm::dd::DdType::Sylvan> const& ddManager, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+
+            template storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> getRowColumnDiagonal(storm::dd::DdManager<storm::dd::DdType::Sylvan> const& ddManager, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+
         }
     }
 }

src/storm/utility/dd.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <set>
+#include <vector>
 
 #include "storm/storage/dd/DdType.h"
 
@@ -9,8 +10,14 @@ namespace storm {
         class Variable;
     }
     namespace dd {
+        template<storm::dd::DdType Type>
+        class DdManager;
+        
         template<storm::dd::DdType Type>
         class Bdd;
+
+        template<storm::dd::DdType Type, typename ValueType>
+        class Add;
     }
     
     namespace utility {
@@ -19,6 +26,12 @@ namespace storm {
             template <storm::dd::DdType Type>
             storm::dd::Bdd<Type> computeReachableStates(storm::dd::Bdd<Type> const& initialStates, storm::dd::Bdd<Type> const& transitions, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables);
             
+            template <storm::dd::DdType Type, typename ValueType>
+            storm::dd::Add<Type, ValueType> getRowColumnDiagonal(storm::dd::DdManager<Type> const& ddManager, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+
+            template <storm::dd::DdType Type>
+            storm::dd::Bdd<Type> getRowColumnDiagonal(storm::dd::DdManager<Type> const& ddManager, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+                        
         }
     }
 }

src/storm/utility/solver.cpp

@@ -3,6 +3,7 @@
 #include <vector>
 
 #include "storm/solver/SymbolicLinearEquationSolver.h"
+#include "storm/solver/SymbolicEliminationLinearEquationSolver.h"
 #include "storm/solver/SymbolicMinMaxLinearEquationSolver.h"
 #include "storm/solver/SymbolicGameSolver.h"
 #include "storm/solver/GameSolver.h"
@@ -27,7 +28,14 @@ namespace storm {
             
             template<storm::dd::DdType Type, typename ValueType>
             std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>> SymbolicLinearEquationSolverFactory<Type, ValueType>::create(storm::dd::Add<Type, ValueType> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
-                return std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>>(new storm::solver::SymbolicLinearEquationSolver<Type, ValueType>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs));
+                
+                storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::CoreSettings>().getEquationSolver();
+                switch (equationSolver) {
+                    case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<Type, ValueType>>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                        break;
+                    default:
+                        return std::make_unique<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                }
             }
             
             template<storm::dd::DdType Type, typename ValueType>
@@ -71,7 +79,7 @@ namespace storm {
             std::unique_ptr<storm::solver::LpSolver> GurobiLpSolverFactory::create(std::string const& name) const {
                 return LpSolverFactory::create(name, storm::solver::LpSolverTypeSelection::Gurobi);
             }
-
+            
             std::unique_ptr<storm::solver::LpSolver> Z3LpSolverFactory::create(std::string const& name) const {
                 return LpSolverFactory::create(name, storm::solver::LpSolverTypeSelection::Z3);
             }