made elimination-based linear solver work in an alpha version. changed minor things in Eigen's SparseLU implementation to make it work with rational numbers and rational functions

Former-commit-id: e5622bd981

resources/3rdparty/eigen-3.2.6/Eigen/src/SparseCore/SparseMatrix.h

@@ -384,7 +384,7 @@ class SparseMatrix
       eigen_assert( (m_outerIndex[outer+1]-m_outerIndex[outer]==0 || m_data.index(m_data.size()-1)<inner) && "Invalid ordered insertion (invalid inner index)");
       Index p = m_outerIndex[outer+1];
       ++m_outerIndex[outer+1];
-      m_data.append(0, inner);
+      m_data.append(Scalar(0), inner);
       return m_data.value(p);
     }
 
@@ -844,7 +844,7 @@ public:
 
       Index p = m_outerIndex[outer] + m_innerNonZeros[outer]++;
       m_data.index(p) = inner;
-      return (m_data.value(p) = 0);
+      return (m_data.value(p) = Scalar(0));
     }
 
 private:

resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU.h

@@ -86,7 +86,7 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ
     typedef internal::SparseLUImpl<Scalar, Index> Base;
     
   public:
-    SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
+    SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1),m_detPermR(1)
     {
       initperfvalues(); 
     }

resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_column_bmod.h

@@ -128,7 +128,7 @@ Index SparseLUImpl<Scalar,Index>::column_bmod(const Index jcol, const Index nseg
   {
     irow = glu.lsub(isub);
     glu.lusup(nextlu) = dense(irow);
-    dense(irow) = Scalar(0.0); 
+    dense(irow) = Scalar(0); 
     ++nextlu; 
   }
   

resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h

@@ -86,7 +86,7 @@ Index SparseLUImpl<Scalar,Index>::copy_to_ucol(const Index jcol, const Index nse
           irow = glu.lsub(isub); 
           glu.usub(nextu) = perm_r(irow); // Unlike the L part, the U part is stored in its final order
           glu.ucol(nextu) = dense(irow); 
-          dense(irow) = Scalar(0.0); 
+          dense(irow) = Scalar(0); 
           nextu++;
           isub++;
         }

resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_panel_bmod.h

@@ -122,7 +122,7 @@ void SparseLUImpl<Scalar,Index>::panel_bmod(const Index m, const Index w, const
         
         Index isub = lptr + no_zeros;
         Index off = u_rows-segsize;
-        for (Index i = 0; i < off; i++) U(i,u_col) = 0;
+        for (Index i = 0; i < off; i++) U(i,u_col) = Scalar(0);
         for (Index i = 0; i < segsize; i++)
         {
           Index irow = glu.lsub(isub); 
@@ -172,7 +172,7 @@ void SparseLUImpl<Scalar,Index>::panel_bmod(const Index m, const Index w, const
         {
           Index irow = glu.lsub(isub++); 
           dense_col(irow) = U.coeff(i+off,u_col);
-          U.coeffRef(i+off,u_col) = 0;
+          U.coeffRef(i+off,u_col) = Scalar(0);
         }
         
         // Scatter l into SPA dense[]
@@ -180,7 +180,7 @@ void SparseLUImpl<Scalar,Index>::panel_bmod(const Index m, const Index w, const
         {
           Index irow = glu.lsub(isub++); 
           dense_col(irow) -= L.coeff(i,u_col);
-          L.coeffRef(i,u_col) = 0;
+          L.coeffRef(i,u_col) = Scalar(0);
         }
         u_col++;
       }

resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_pivotL.h

@@ -71,42 +71,43 @@ Index SparseLUImpl<Scalar,Index>::pivotL(const Index jcol, const RealScalar& dia
   
   // Determine the largest abs numerical value for partial pivoting 
   Index diagind = iperm_c(jcol); // diagonal index 
-  RealScalar pivmax(-1.0);
-  Index pivptr = nsupc; 
+//  RealScalar pivmax(-1);
+  Index pivptr = nsupc;
   Index diag = emptyIdxLU; 
-  RealScalar rtemp;
-  Index isub, icol, itemp, k; 
+//  RealScalar rtemp;
+  Index isub, icol, itemp, k;
   for (isub = nsupc; isub < nsupr; ++isub) {
-    using std::abs;
-    rtemp = abs(lu_col_ptr[isub]);
-    if (rtemp > pivmax) {
-      pivmax = rtemp; 
+//    using std::abs;
+//    rtemp = abs(lu_col_ptr[isub]);
+//    if (rtemp > pivmax) {
+//      pivmax = rtemp;
       pivptr = isub;
-    } 
+//    } 
     if (lsub_ptr[isub] == diagind) diag = isub;
   }
   
   // Test for singularity
-  if ( pivmax <= RealScalar(0.0) ) {
-    // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero
-    pivrow = pivmax < RealScalar(0.0) ? diagind : lsub_ptr[pivptr];
-    perm_r(pivrow) = jcol;
-    return (jcol+1);
-  }
+//  if ( pivmax <= RealScalar(0) ) {
+//    // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero
+//    pivrow = pivmax < RealScalar(0) ? diagind : lsub_ptr[pivptr];
+//    perm_r(pivrow) = jcol;
+//    return (jcol+1);
+//  }
   
-  RealScalar thresh = diagpivotthresh * pivmax; 
+//  RealScalar thresh = diagpivotthresh * pivmax; 
   
   // Choose appropriate pivotal element 
   
   {
     // Test if the diagonal element can be used as a pivot (given the threshold value)
-    if (diag >= 0 ) 
-    {
+//    if (diag >= 0 ) 
+//    {
       // Diagonal element exists
-      using std::abs;
-      rtemp = abs(lu_col_ptr[diag]);
-      if (rtemp != 0.0 && rtemp >= thresh) pivptr = diag;
-    }
+//      using std::abs;
+//      rtemp = abs(lu_col_ptr[diag]);
+//      if (rtemp != 0.0 && rtemp >= thresh)
+          pivptr = diag;
+//    }
     pivrow = lsub_ptr[pivptr];
   }
   
@@ -125,7 +126,7 @@ Index SparseLUImpl<Scalar,Index>::pivotL(const Index jcol, const RealScalar& dia
     }
   }
   // cdiv operations
-  Scalar temp = Scalar(1.0) / lu_col_ptr[nsupc];
+  Scalar temp = Scalar(1) / lu_col_ptr[nsupc];
   for (k = nsupc+1; k < nsupr; k++)
     lu_col_ptr[k] *= temp; 
   return 0;

src/adapters/EigenAdapter.cpp

@@ -21,5 +21,7 @@ namespace storm {
         }
         
         template std::unique_ptr<Eigen::SparseMatrix<double>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<double> const& matrix);
+        template std::unique_ptr<Eigen::SparseMatrix<storm::RationalNumber>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<storm::RationalNumber> const& matrix);
+        template std::unique_ptr<Eigen::SparseMatrix<storm::RationalFunction>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix);
     }
 }

src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp

@@ -19,7 +19,7 @@
 #include "src/logic/FragmentSpecification.h"
 
 #include "src/solver/stateelimination/LongRunAverageEliminator.h"
-#include "src/solver/stateelimination/ConditionalEliminator.h"
+#include "src/solver/stateelimination/ConditionalStateEliminator.h"
 #include "src/solver/stateelimination/PrioritizedStateEliminator.h"
 #include "src/solver/stateelimination/StaticStatePriorityQueue.h"
 #include "src/solver/stateelimination/DynamicStatePriorityQueue.h"
@@ -685,7 +685,7 @@ namespace storm {
             STORM_LOG_INFO("Eliminating " << numberOfStatesToEliminate << " states using the state elimination technique." << std::endl);
             performPrioritizedStateElimination(statePriorities, flexibleMatrix, flexibleBackwardTransitions, oneStepProbabilities, this->getModel().getInitialStates(), true);
             
-            storm::solver::stateelimination::ConditionalEliminator<ValueType> stateEliminator = storm::solver::stateelimination::ConditionalEliminator<ValueType>(flexibleMatrix, flexibleBackwardTransitions, oneStepProbabilities, phiStates, psiStates);
+            storm::solver::stateelimination::ConditionalStateEliminator<ValueType> stateEliminator = storm::solver::stateelimination::ConditionalStateEliminator<ValueType>(flexibleMatrix, flexibleBackwardTransitions, oneStepProbabilities, phiStates, psiStates);
             
             // Eliminate the transitions going into the initial state (if there are any).
             if (!flexibleBackwardTransitions.getRow(*newInitialStates.begin()).empty()) {
@@ -948,7 +948,7 @@ namespace storm {
             // Finally, eliminate the entry states (if we are required to do so).
             if (eliminateEntryStates) {
                 STORM_LOG_TRACE("Finally, eliminating entry states.");
-                std::shared_ptr<StatePriorityQueue> naivePriorities = createNaivePriorityQueue(entryStates);
+                std::shared_ptr<StatePriorityQueue> naivePriorities = createStatePriorityQueue(entryStates);
                 performPrioritizedStateElimination(naivePriorities, matrix, backwardTransitions, values, initialStates, computeResultsForInitialStatesOnly);
                 STORM_LOG_TRACE("Eliminated/added entry states.");
             } else {

src/settings/modules/MarkovChainSettings.cpp

@@ -41,9 +41,9 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, engineOptionName, false, "Sets which engine is used for model building and model checking.").setShortName(engineOptionShortName)
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the engine to use. Available are {sparse, hybrid, dd, expl, abs}.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(engines)).setDefaultValueString("sparse").build()).build());
                 
-                std::vector<std::string> linearEquationSolver = {"gmm++", "native", "eigen"};
+                std::vector<std::string> linearEquationSolver = {"gmm++", "native", "eigen", "elimination"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, eqSolverOptionName, false, "Sets which solver is preferred for solving systems of linear equations.")
-                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the solver to prefer. Available are: gmm++, native and eigen.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(linearEquationSolver)).setDefaultValueString("gmm++").build()).build());
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the solver to prefer. Available are: gmm++, native, eigen, elimination.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(linearEquationSolver)).setDefaultValueString("gmm++").build()).build());
                 
                 std::vector<std::string> ddLibraries = {"cudd", "sylvan"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, ddLibraryOptionName, false, "Sets which library is preferred for decision-diagram operations.")
@@ -87,6 +87,8 @@ namespace storm {
                     return storm::solver::EquationSolverType::Native;
                 } else if (equationSolverName == "eigen") {
                     return storm::solver::EquationSolverType::Eigen;
+                } else if (equationSolverName == "elimination") {
+                    return storm::solver::EquationSolverType::Elimination;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown equation solver '" << equationSolverName << "'.");
             }

src/solver/EigenLinearEquationSolver.cpp

@@ -9,7 +9,7 @@ namespace storm {
     namespace solver {
      
         template<typename ValueType>
-        EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A)), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner) {
+        EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A)), method(method), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner) {
             // Intentionally left empty.
         }
         
@@ -44,31 +44,31 @@ namespace storm {
         template<typename ValueType>
         void EigenLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             // Translate the vectors x and b into Eigen's format.
-            Eigen::VectorXd eigenB(b.size());
+            Eigen::Matrix<ValueType, Eigen::Dynamic, 1> eigenB(b.size());
             for (uint64_t index = 0; index < b.size(); ++index) {
                 eigenB[index] = b[index];
             }
             
-            Eigen::VectorXd eigenX(x.size());
+            Eigen::Matrix<ValueType, Eigen::Dynamic, 1> eigenX(x.size());
             for (uint64_t index = 0; index < x.size(); ++index) {
                 eigenX[index] = x[index];
             }
             
             if (method == SolutionMethod::SparseLU) {
-                Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
+                Eigen::SparseLU<Eigen::SparseMatrix<ValueType>, Eigen::COLAMDOrdering<int>> solver;
                 solver.compute(*eigenA);
                 solver._solve(eigenB, eigenX);
             } else {
                 if (preconditioner == Preconditioner::Ilu) {
-                    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::IncompleteLUT<double>> solver;
+                    Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::IncompleteLUT<ValueType>> solver;
                     solver.compute(*eigenA);
                     solver._solve(eigenB, eigenX);
                 } else if (preconditioner == Preconditioner::Diagonal) {
-                    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::DiagonalPreconditioner<double>> solver;
+                    Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::DiagonalPreconditioner<ValueType>> solver;
                     solver.compute(*eigenA);
                     solver._solve(eigenB, eigenX);
                 } else {
-                    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::IdentityPreconditioner> solver;
+                    Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::IdentityPreconditioner> solver;
                     solver.compute(*eigenA);
                     solver._solve(eigenB, eigenX);
                 }
@@ -82,16 +82,73 @@ namespace storm {
         
         template<typename ValueType>
         void EigenLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+            // Translate the vectors x and b into Eigen's format.
+            std::unique_ptr<Eigen::Matrix<ValueType, Eigen::Dynamic, 1>> eigenB;
+            if (b != nullptr) {
+                eigenB = std::make_unique<Eigen::Matrix<ValueType, Eigen::Dynamic, 1>>(b->size());
+                for (uint64_t index = 0; index < b->size(); ++index) {
+                    (*eigenB)[index] = (*b)[index];
+                }
+            }
+            Eigen::Matrix<ValueType, Eigen::Dynamic, 1> eigenX(x.size());
+            for (uint64_t index = 0; index < x.size(); ++index) {
+                eigenX[index] = x[index];
+            }
+            
+            // Perform n matrix-vector multiplications.
+            for (uint64_t iteration = 0; iteration < n; ++iteration) {
+                eigenX = *eigenA * eigenX;
+                if (eigenB != nullptr) {
+                    eigenX += *eigenB;
+                }
+            }
+            
+            // Translate the solution from Eigen's format into our representation.
+            for (uint64_t index = 0; index < eigenX.size(); ++index) {
+                x[index] = eigenX[index];
+            }
+        }
+        
+        // Specialization form storm::RationalNumber
+        
+        EigenLinearEquationSolver<storm::RationalNumber>::EigenLinearEquationSolver(storm::storage::SparseMatrix<storm::RationalNumber> const& A, SolutionMethod method) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<storm::RationalNumber>(A)), method(method) {
+            // Intentionally left empty.
+        }
+        
+        void EigenLinearEquationSolver<storm::RationalNumber>::solveEquationSystem(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b, std::vector<storm::RationalNumber>* multiplyResult) const {
+            // Translate the vectors x and b into Eigen's format.
+            Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> eigenB(b.size());
+            for (uint64_t index = 0; index < b.size(); ++index) {
+                eigenB[index] = b[index];
+            }
+            
+            Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> eigenX(x.size());
+            for (uint64_t index = 0; index < x.size(); ++index) {
+                eigenX[index] = x[index];
+            }
             
+            if (method == SolutionMethod::SparseLU) {
+                Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalNumber>, Eigen::COLAMDOrdering<int>> solver;
+                solver.compute(*eigenA);
+                solver._solve(eigenB, eigenX);
+            }
+            
+            // Translate the solution from Eigen's format into our representation.
+            for (uint64_t index = 0; index < eigenX.size(); ++index) {
+                x[index] = eigenX[index];
+            }
+        }
+        
+        void EigenLinearEquationSolver<storm::RationalNumber>::performMatrixVectorMultiplication(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const* b, uint_fast64_t n, std::vector<storm::RationalNumber>* multiplyResult) const {
             // Translate the vectors x and b into Eigen's format.
-            std::unique_ptr<Eigen::VectorXd> eigenB;
+            std::unique_ptr<Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>> eigenB;
             if (b != nullptr) {
-                eigenB = std::make_unique<Eigen::VectorXd>(b->size());
+                eigenB = std::make_unique<Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>>(b->size());
                 for (uint64_t index = 0; index < b->size(); ++index) {
                     (*eigenB)[index] = (*b)[index];
                 }
             }
-            Eigen::VectorXd eigenX(x.size());
+            Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> eigenX(x.size());
             for (uint64_t index = 0; index < x.size(); ++index) {
                 eigenX[index] = x[index];
             }
@@ -110,8 +167,67 @@ namespace storm {
             }
         }
         
-        template class EigenLinearEquationSolver<double>;
+        // Specialization form storm::RationalFunction
         
+        EigenLinearEquationSolver<storm::RationalFunction>::EigenLinearEquationSolver(storm::storage::SparseMatrix<storm::RationalFunction> const& A, SolutionMethod method) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<storm::RationalFunction>(A)), method(method) {
+            // Intentionally left empty.
+        }
+        
+        void EigenLinearEquationSolver<storm::RationalFunction>::solveEquationSystem(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b, std::vector<storm::RationalFunction>* multiplyResult) const {
+            // Translate the vectors x and b into Eigen's format.
+            Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> eigenB(b.size());
+            for (uint64_t index = 0; index < b.size(); ++index) {
+                eigenB[index] = b[index];
+            }
+            
+            Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> eigenX(x.size());
+            for (uint64_t index = 0; index < x.size(); ++index) {
+                eigenX[index] = x[index];
+            }
+            
+            if (method == SolutionMethod::SparseLU) {
+                Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalFunction>, Eigen::COLAMDOrdering<int>> solver;
+                solver.compute(*eigenA);
+                solver._solve(eigenB, eigenX);
+            }
+            
+            // Translate the solution from Eigen's format into our representation.
+            for (uint64_t index = 0; index < eigenX.size(); ++index) {
+                x[index] = eigenX[index];
+            }
+        }
+        
+        void EigenLinearEquationSolver<storm::RationalFunction>::performMatrixVectorMultiplication(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const* b, uint_fast64_t n, std::vector<storm::RationalFunction>* multiplyResult) const {
+            // Translate the vectors x and b into Eigen's format.
+            std::unique_ptr<Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>> eigenB;
+            if (b != nullptr) {
+                eigenB = std::make_unique<Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>>(b->size());
+                for (uint64_t index = 0; index < b->size(); ++index) {
+                    (*eigenB)[index] = (*b)[index];
+                }
+            }
+            Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> eigenX(x.size());
+            for (uint64_t index = 0; index < x.size(); ++index) {
+                eigenX[index] = x[index];
+            }
+            
+            // Perform n matrix-vector multiplications.
+            for (uint64_t iteration = 0; iteration < n; ++iteration) {
+                eigenX = *eigenA * eigenX;
+                if (eigenB != nullptr) {
+                    eigenX += *eigenB;
+                }
+            }
+            
+            // Translate the solution from Eigen's format into our representation.
+            for (uint64_t index = 0; index < eigenX.size(); ++index) {
+                x[index] = eigenX[index];
+            }
+        }
+        
+        template class EigenLinearEquationSolver<double>;
+        template class EigenLinearEquationSolver<storm::RationalNumber>;
+//        template class EigenLinearEquationSolver<storm::RationalFunction>;
         
     }
 }

src/solver/EigenLinearEquationSolver.h

@@ -34,7 +34,7 @@ namespace storm {
             storm::storage::SparseMatrix<ValueType> const* originalA;
             
             // The (eigen) matrix associated with this equation solver.
-            std::unique_ptr<Eigen::SparseMatrix<double>> eigenA;
+            std::unique_ptr<Eigen::SparseMatrix<ValueType>> eigenA;
 
             // The method to use for solving linear equation systems.
             SolutionMethod method;
@@ -49,5 +49,52 @@ namespace storm {
             Preconditioner preconditioner;
         };
         
+        template<>
+        class EigenLinearEquationSolver<storm::RationalNumber> : public LinearEquationSolver<storm::RationalNumber> {
+        public:
+            enum class SolutionMethod {
+                SparseLU
+            };
+            
+            EigenLinearEquationSolver(storm::storage::SparseMatrix<storm::RationalNumber> const& A, SolutionMethod method = SolutionMethod::SparseLU);
+            
+            virtual void solveEquationSystem(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b, std::vector<storm::RationalNumber>* multiplyResult = nullptr) const override;
+            
+            virtual void performMatrixVectorMultiplication(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const* b, uint_fast64_t n = 1, std::vector<storm::RationalNumber>* multiplyResult = nullptr) const override;
+            
+        private:
+            // A pointer to the original sparse matrix given to this solver.
+            storm::storage::SparseMatrix<storm::RationalNumber> const* originalA;
+            
+            // The (eigen) matrix associated with this equation solver.
+            std::unique_ptr<Eigen::SparseMatrix<storm::RationalNumber>> eigenA;
+            
+            // The method to use for solving linear equation systems.
+            SolutionMethod method;
+        };
+        
+        template<>
+        class EigenLinearEquationSolver<storm::RationalFunction> : public LinearEquationSolver<storm::RationalFunction> {
+        public:
+            enum class SolutionMethod {
+                SparseLU
+            };
+            
+            EigenLinearEquationSolver(storm::storage::SparseMatrix<storm::RationalFunction> const& A, SolutionMethod method = SolutionMethod::SparseLU);
+            
+            virtual void solveEquationSystem(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b, std::vector<storm::RationalFunction>* multiplyResult = nullptr) const override;
+            
+            virtual void performMatrixVectorMultiplication(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const* b, uint_fast64_t n = 1, std::vector<storm::RationalFunction>* multiplyResult = nullptr) const override;
+            
+        private:
+            // A pointer to the original sparse matrix given to this solver.
+            storm::storage::SparseMatrix<storm::RationalFunction> const* originalA;
+            
+            // The (eigen) matrix associated with this equation solver.
+            std::unique_ptr<Eigen::SparseMatrix<storm::RationalFunction>> eigenA;
+            
+            // The method to use for solving linear equation systems.
+            SolutionMethod method;
+        };
     }
 }

src/solver/EliminationLinearEquationSolver.cpp

@@ -1,9 +1,19 @@
 #include "src/solver/EliminationLinearEquationSolver.h"
 
+#include <numeric>
+
+#include "src/solver/stateelimination/StatePriorityQueue.h"
+#include "src/solver/stateelimination/PrioritizedStateEliminator.h"
+
+#include "src/utility/macros.h"
 #include "src/utility/vector.h"
+#include "src/utility/stateelimination.h"
 
 namespace storm {
     namespace solver {
+        
+        using namespace stateelimination;
+        
         template<typename ValueType>
         EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : A(A) {
             // Intentionally left empty.
@@ -11,7 +21,54 @@ namespace storm {
         
         template<typename ValueType>
         void EliminationLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
-            // TODO: implement state-elimination here.
+            STORM_LOG_WARN_COND(multiplyResult == nullptr, "Providing scratch memory will not improve the performance of this solver.");
+            
+//            std::cout << "input:" << std::endl;
+//            std::cout << "A:" << std::endl;
+//            std::cout << A << std::endl;
+//            std::cout << "b:" << std::endl;
+//            for (auto const& e : b) {
+//                std::cout << e << std::endl;
+//            }
+            
+            // Create a naive priority queue.
+            // TODO: improve.
+            std::vector<storm::storage::sparse::state_type> allRows(x.size());
+            std::iota(allRows.begin(), allRows.end(), 0);
+            std::shared_ptr<StatePriorityQueue> priorityQueue = storm::utility::stateelimination::createStatePriorityQueue(allRows);
+            
+            // Initialize the solution to the right-hand side of the equation system.
+            x = b;
+            
+            // Translate the matrix and its transpose into the flexible format.
+            // We need to revert the matrix from the equation system format to the probability matrix format. That is,
+            // from (I-P), we construct P. Note that for the backwards transitions, this does not need to be done, as all
+            // entries are set to one anyway.
+            storm::storage::FlexibleSparseMatrix<ValueType> flexibleMatrix(A, false, true);
+            storm::storage::FlexibleSparseMatrix<ValueType> flexibleBackwardTransitions(A.transpose(), true, true);
+            
+//            std::cout << "intermediate:" << std::endl;
+//            std::cout << "flexibleMatrix:" << std::endl;
+//            std::cout << flexibleMatrix << std::endl;
+//            std::cout << "backward:" << std::endl;
+//            std::cout << flexibleBackwardTransitions << std::endl;
+            
+            // Create a state eliminator to perform the actual elimination.
+            PrioritizedStateEliminator<ValueType> eliminator(flexibleMatrix, flexibleBackwardTransitions, priorityQueue, x);
+            
+            std::cout << "eliminating" << std::endl;
+            while (priorityQueue->hasNext()) {
+                auto state = priorityQueue->pop();
+                eliminator.eliminateState(state, false);
+            }
+            
+//            std::cout << "output:" << std::endl;
+//            std::cout << "x:" << std::endl;
+//            for (auto const& e : x) {
+//                std::cout << e << std::endl;
+//            }
+            
+            std::cout << "done" << std::endl;
         }
         
         template<typename ValueType>
@@ -52,8 +109,6 @@ namespace storm {
         }
 
         template class EliminationLinearEquationSolver<double>;
-        
-        // TODO: make this work with the proper implementation of solveEquationSystem.
         template class EliminationLinearEquationSolver<storm::RationalNumber>;
         template class EliminationLinearEquationSolver<storm::RationalFunction>;
         

src/solver/SolverSelectionOptions.h

@@ -9,7 +9,7 @@ namespace storm {
         ExtendEnumsWithSelectionField(MinMaxTechnique, PolicyIteration, ValueIteration)
         
         ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk)
-        ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Topological)
+        ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination, Topological)
         ExtendEnumsWithSelectionField(SmtSolverType, Z3, Mathsat)
     }
 } 

src/solver/stateelimination/ConditionalEliminator.cpp → src/solver/stateelimination/ConditionalStateEliminator.cpp

@@ -1,4 +1,4 @@
-#include "src/solver/stateelimination/ConditionalEliminator.h"
+#include "src/solver/stateelimination/ConditionalStateEliminator.h"
 
 #include "src/utility/macros.h"
 #include "src/utility/constants.h"
@@ -8,26 +8,21 @@ namespace storm {
         namespace stateelimination {
             
             template<typename ValueType>
-            ConditionalEliminator<ValueType>::ConditionalEliminator(storm::storage::FlexibleSparseMatrix<ValueType>& transitionMatrix, storm::storage::FlexibleSparseMatrix<ValueType>& backwardTransitions, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector& phiStates, storm::storage::BitVector& psiStates) : StateEliminator<ValueType>(transitionMatrix, backwardTransitions), oneStepProbabilities(oneStepProbabilities), phiStates(phiStates), psiStates(psiStates), filterLabel(StateLabel::NONE) {
+            ConditionalStateEliminator<ValueType>::ConditionalStateEliminator(storm::storage::FlexibleSparseMatrix<ValueType>& transitionMatrix, storm::storage::FlexibleSparseMatrix<ValueType>& backwardTransitions, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector& phiStates, storm::storage::BitVector& psiStates) : StateEliminator<ValueType>(transitionMatrix, backwardTransitions), oneStepProbabilities(oneStepProbabilities), phiStates(phiStates), psiStates(psiStates), filterLabel(StateLabel::NONE) {
             }
             
             template<typename ValueType>
-            void ConditionalEliminator<ValueType>::updateValue(storm::storage::sparse::state_type const& state, ValueType const& loopProbability) {
+            void ConditionalStateEliminator<ValueType>::updateValue(storm::storage::sparse::state_type const& state, ValueType const& loopProbability) {
                 oneStepProbabilities[state] = storm::utility::simplify(loopProbability * oneStepProbabilities[state]);
             }
             
             template<typename ValueType>
-            void ConditionalEliminator<ValueType>::updatePredecessor(storm::storage::sparse::state_type const& predecessor, ValueType const& probability, storm::storage::sparse::state_type const& state) {
+            void ConditionalStateEliminator<ValueType>::updatePredecessor(storm::storage::sparse::state_type const& predecessor, ValueType const& probability, storm::storage::sparse::state_type const& state) {
                 oneStepProbabilities[predecessor] = storm::utility::simplify(oneStepProbabilities[predecessor] * storm::utility::simplify(probability * oneStepProbabilities[state]));
             }
-            
-            template<typename ValueType>
-            void ConditionalEliminator<ValueType>::updatePriority(storm::storage::sparse::state_type const& state) {
-                // Do nothing
-            }
-            
+                        
             template<typename ValueType>
-            bool ConditionalEliminator<ValueType>::filterPredecessor(storm::storage::sparse::state_type const& state) {
+            bool ConditionalStateEliminator<ValueType>::filterPredecessor(storm::storage::sparse::state_type const& state) {
                 // TODO find better solution than flag
                 switch (filterLabel) {
                     case StateLabel::PHI:
@@ -41,34 +36,34 @@ namespace storm {
             }
             
             template<typename ValueType>
-            bool ConditionalEliminator<ValueType>::isFilterPredecessor() const {
+            bool ConditionalStateEliminator<ValueType>::isFilterPredecessor() const {
                 return true;
             }
             
             template<typename ValueType>
-            void ConditionalEliminator<ValueType>::setFilterPhi() {
+            void ConditionalStateEliminator<ValueType>::setFilterPhi() {
                 filterLabel = StateLabel::PHI;
             }
             
             template<typename ValueType>
-            void ConditionalEliminator<ValueType>::setFilterPsi() {
+            void ConditionalStateEliminator<ValueType>::setFilterPsi() {
                 filterLabel = StateLabel::PSI;
             }
             
             template<typename ValueType>
-            void ConditionalEliminator<ValueType>::setFilter(StateLabel const& stateLabel) {
+            void ConditionalStateEliminator<ValueType>::setFilter(StateLabel const& stateLabel) {
                 filterLabel = stateLabel;
             }
             
             template<typename ValueType>
-            void ConditionalEliminator<ValueType>::unsetFilter() {
+            void ConditionalStateEliminator<ValueType>::unsetFilter() {
                 filterLabel = StateLabel::NONE;
             }
             
-            template class ConditionalEliminator<double>;
+            template class ConditionalStateEliminator<double>;
             
 #ifdef STORM_HAVE_CARL
-            template class ConditionalEliminator<storm::RationalFunction>;
+            template class ConditionalStateEliminator<storm::RationalFunction>;
 #endif
             
         } // namespace stateelimination

src/solver/stateelimination/ConditionalEliminator.h → src/solver/stateelimination/ConditionalStateEliminator.h

@@ -1,5 +1,5 @@
-#ifndef STORM_SOLVER_STATEELIMINATION_CONDITIONALELIMINATOR_H_
-#define STORM_SOLVER_STATEELIMINATION_CONDITIONALELIMINATOR_H_
+#ifndef STORM_SOLVER_STATEELIMINATION_CONDITIONALSTATEELIMINATOR_H_
+#define STORM_SOLVER_STATEELIMINATION_CONDITIONALSTATEELIMINATOR_H_
 
 #include "src/solver/stateelimination/StateEliminator.h"
 
@@ -10,16 +10,15 @@ namespace storm {
         namespace stateelimination {
             
             template<typename ValueType>
-            class ConditionalEliminator : public StateEliminator<ValueType> {
+            class ConditionalStateEliminator : public StateEliminator<ValueType> {
             public:
                 enum class StateLabel { NONE, PHI, PSI };
                 
-                ConditionalEliminator(storm::storage::FlexibleSparseMatrix<ValueType>& transitionMatrix, storm::storage::FlexibleSparseMatrix<ValueType>& backwardTransitions, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector& phiStates, storm::storage::BitVector& psiStates);
+                ConditionalStateEliminator(storm::storage::FlexibleSparseMatrix<ValueType>& transitionMatrix, storm::storage::FlexibleSparseMatrix<ValueType>& backwardTransitions, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector& phiStates, storm::storage::BitVector& psiStates);
                 
                 // Instantiaton of Virtual methods
                 void updateValue(storm::storage::sparse::state_type const& state, ValueType const& loopProbability) override;
                 void updatePredecessor(storm::storage::sparse::state_type const& predecessor, ValueType const& probability, storm::storage::sparse::state_type const& state) override;
-                void updatePriority(storm::storage::sparse::state_type const& state) override;
                 bool filterPredecessor(storm::storage::sparse::state_type const& state) override;
                 bool isFilterPredecessor() const override;
                 
@@ -39,4 +38,4 @@ namespace storm {
     } // namespace storage
 } // namespace storm
 
-#endif // STORM_SOLVER_STATEELIMINATION_CONDITIONALELIMINATOR_H_
+#endif // STORM_SOLVER_STATEELIMINATION_CONDITIONALSTATEELIMINATOR_H_

src/solver/stateelimination/PrioritizedStateEliminator.cpp

@@ -29,7 +29,8 @@ namespace storm {
             }
             
             template class PrioritizedStateEliminator<double>;
-            
+            template class PrioritizedStateEliminator<storm::RationalNumber>;
+
 #ifdef STORM_HAVE_CARL
             template class PrioritizedStateEliminator<storm::RationalFunction>;
 #endif

src/solver/stateelimination/StateEliminator.cpp

@@ -264,6 +264,7 @@ namespace storm {
             }
             
             template class StateEliminator<double>;
+            template class StateEliminator<storm::RationalNumber>;
             
 #ifdef STORM_HAVE_CARL
             template class StateEliminator<storm::RationalFunction>;

src/storage/FlexibleSparseMatrix.cpp

@@ -3,7 +3,10 @@
 #include "src/storage/SparseMatrix.h"
 #include "src/storage/BitVector.h"
 #include "src/adapters/CarlAdapter.h"
+
+#include "src/utility/macros.h"
 #include "src/utility/constants.h"
+#include "src/exceptions/InvalidArgumentException.h"
 
 namespace storm {
     namespace storage {
@@ -13,7 +16,9 @@ namespace storm {
         }
         
         template<typename ValueType>
-        FlexibleSparseMatrix<ValueType>::FlexibleSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix, bool setAllValuesToOne) : data(matrix.getRowCount()), columnCount(matrix.getColumnCount()), nonzeroEntryCount(matrix.getNonzeroEntryCount()), trivialRowGrouping(matrix.hasTrivialRowGrouping()) {
+        FlexibleSparseMatrix<ValueType>::FlexibleSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix, bool setAllValuesToOne, bool revertEquationSystem) : data(matrix.getRowCount()), columnCount(matrix.getColumnCount()), nonzeroEntryCount(matrix.getNonzeroEntryCount()), trivialRowGrouping(matrix.hasTrivialRowGrouping()) {
+            STORM_LOG_THROW(!revertEquationSystem || trivialRowGrouping, storm::exceptions::InvalidArgumentException, "Illegal option for creating flexible matrix.");
+            
             if (!trivialRowGrouping) {
                 rowGroupIndices = matrix.getRowGroupIndices();
             }
@@ -23,12 +28,31 @@ namespace storm {
                 for (auto const& element : row) {
                     // If the probability is zero, we skip this entry.
                     if (storm::utility::isZero(element.getValue())) {
-                        continue;
+                        if (revertEquationSystem && rowIndex == element.getColumn()) {
+                            getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one<ValueType>());
+                        } else {
+                            continue;
+                        }
                     }
                     if (setAllValuesToOne) {
-                        getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one<ValueType>());
+                        if (revertEquationSystem && element.getColumn() == rowIndex && storm::utility::isOne(element.getValue())) {
+                            continue;
+                        } else {
+                            getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one<ValueType>());
+                        }
                     } else {
-                        getRow(rowIndex).emplace_back(element);
+                        if (revertEquationSystem) {
+                            if (element.getColumn() == rowIndex) {
+                                if (storm::utility::isOne(element.getValue())) {
+                                    continue;
+                                }
+                                getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one<ValueType>() - element.getValue());
+                            } else {
+                                getRow(rowIndex).emplace_back(element.getColumn(), -element.getValue());
+                            }
+                        } else {
+                            getRow(rowIndex).emplace_back(element);
+                        }
                     }
                 }
             }
@@ -245,7 +269,9 @@ namespace storm {
         // Explicitly instantiate the matrix.
         template class FlexibleSparseMatrix<double>;
         template std::ostream& operator<<(std::ostream& out, FlexibleSparseMatrix<double> const& matrix);
-        
+        template class FlexibleSparseMatrix<storm::RationalNumber>;
+        template std::ostream& operator<<(std::ostream& out, FlexibleSparseMatrix<storm::RationalNumber> const& matrix);
+
 #ifdef STORM_HAVE_CARL
         template class FlexibleSparseMatrix<storm::RationalFunction>;
         template std::ostream& operator<<(std::ostream& out, FlexibleSparseMatrix<storm::RationalFunction> const& matrix);

src/storage/FlexibleSparseMatrix.h

@@ -44,8 +44,10 @@ namespace storm {
              * Constructs a flexible sparse matrix from a sparse matrix.
              * @param matrix Sparse matrix to construct from.
              * @param setAllValuesToOne If true, all set entries are set to one. Default is false.
+             * @param revertEquationSystem If true, the matrix that will be created is the matrix (1-A), where A is the
+             * provided matrix.
              */
-            FlexibleSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix, bool setAllValuesToOne = false);
+            FlexibleSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix, bool setAllValuesToOne = false, bool revertEquationSystem = false);
 
             /*!
              * Reserves space for elements in row.

src/utility/constants.cpp

@@ -231,6 +231,24 @@ namespace storm {
         template storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::storage::sparse::state_type>& simplify(storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::storage::sparse::state_type>& matrixEntry);
         template storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::storage::sparse::state_type>&& simplify(storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::storage::sparse::state_type>&& matrixEntry);
 
+        // Instantiations for rational number.
+        template bool isOne(storm::RationalNumber const& value);
+        template bool isZero(storm::RationalNumber const& value);
+        template bool isConstant(storm::RationalNumber const& value);
+        
+        template storm::RationalNumber one();
+        template storm::RationalNumber zero();
+        
+        template double convertNumber(storm::RationalNumber const& number);
+        template storm::RationalNumber convertNumber(double const& number);
+
+//        template storm::RationalNumber pow(storm::RationalNumber const& value, uint_fast64_t exponent);
+        
+        template storm::RationalNumber simplify(storm::RationalNumber value);
+        template storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::RationalNumber> simplify(storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::RationalNumber> matrixEntry);
+        template storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::RationalNumber>& simplify(storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::RationalNumber>& matrixEntry);
+        template storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::RationalNumber>&& simplify(storm::storage::MatrixEntry<storm::storage::sparse::state_type, storm::RationalNumber>&& matrixEntry);
+        
 #ifdef STORM_HAVE_CARL
         template bool isOne(RationalFunction const& value);
         template bool isZero(RationalFunction const& value);
@@ -248,16 +266,6 @@ namespace storm {
         template Polynomial one();
         template Polynomial zero();
         
-        template bool isOne(RationalNumber const& value);
-        template bool isZero(RationalNumber const& value);
-        template bool isConstant(RationalNumber const& value);
-        
-        template RationalNumber one();
-        template RationalNumber zero();
-       
-        template double convertNumber(RationalNumber const& number);
-        template RationalNumber convertNumber(double const& number);
-        
         template bool isOne(Interval const& value);
         template bool isZero(Interval const& value);
         template bool isConstant(Interval const& value);

src/utility/solver.cpp

@@ -8,6 +8,7 @@
 #include "src/solver/NativeLinearEquationSolver.h"
 #include "src/solver/GmmxxLinearEquationSolver.h"
 #include "src/solver/EigenLinearEquationSolver.h"
+#include "src/solver/EliminationLinearEquationSolver.h"
 #include "src/solver/GameSolver.h"
 
 #include "src/solver/NativeMinMaxLinearEquationSolver.h"
@@ -49,13 +50,18 @@ namespace storm {
             std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
                 storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::MarkovChainSettings>().getEquationSolver();
                 switch (equationSolver) {
-                    case storm::solver::EquationSolverType::Gmmxx: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
-                    case storm::solver::EquationSolverType::Native: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix));
-                    case storm::solver::EquationSolverType::Eigen: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::EigenLinearEquationSolver<ValueType>(matrix));
-                    default: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
+                    case storm::solver::EquationSolverType::Gmmxx: return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(matrix);
+                    case storm::solver::EquationSolverType::Native: return std::make_unique<storm::solver::NativeLinearEquationSolver<ValueType>>(matrix);
+                    case storm::solver::EquationSolverType::Eigen: return std::make_unique<storm::solver::EigenLinearEquationSolver<ValueType>>(matrix);
+                    case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::EliminationLinearEquationSolver<ValueType>>(matrix);
+                    default: return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(matrix);
                 }
             }
             
+            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> LinearEquationSolverFactory<storm::RationalFunction>::create(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix) const {
+                return std::make_unique<storm::solver::EigenLinearEquationSolver<storm::RationalFunction>>(matrix);
+            }
+            
             template<typename ValueType>
             std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
                 return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
@@ -205,6 +211,7 @@ namespace storm {
             template class SymbolicGameSolverFactory<storm::dd::DdType::Sylvan, double>;
             template class LinearEquationSolverFactory<double>;
             template class GmmxxLinearEquationSolverFactory<double>;
+            template class EigenLinearEquationSolverFactory<double>;
             template class NativeLinearEquationSolverFactory<double>;
             template class MinMaxLinearEquationSolverFactory<double>;
             template class GameSolverFactory<double>;

src/utility/solver.h

@@ -4,8 +4,10 @@
 #include <set>
 #include <vector>
 #include <memory>
-#include <src/storage/sparse/StateType.h>
 
+#include "src/adapters/CarlAdapter.h"
+
+#include "src/storage/sparse/StateType.h"
 #include "src/storage/dd/DdType.h"
 #include "src/solver/SolverSelectionOptions.h"
 
@@ -87,6 +89,18 @@ namespace storm {
                 virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const;
             };
             
+            template<>
+            class LinearEquationSolverFactory<storm::RationalFunction> {
+            public:
+                /*!
+                 * Creates a new linear equation solver instance with the given matrix.
+                 *
+                 * @param matrix The matrix that defines the equation system.
+                 * @return A pointer to the newly created solver.
+                 */
+                virtual std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> create(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix) const;
+            };
+            
             template<typename ValueType>
             class NativeLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
             public:

src/utility/stateelimination.cpp

@@ -138,16 +138,21 @@ namespace storm {
                 STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Illegal elimination order selected.");
             }
             
-            std::shared_ptr<StatePriorityQueue> createNaivePriorityQueue(storm::storage::BitVector const& states) {
+            std::shared_ptr<StatePriorityQueue> createStatePriorityQueue(storm::storage::BitVector const& states) {
                 std::vector<storm::storage::sparse::state_type> sortedStates(states.begin(), states.end());
                 return std::make_shared<StaticStatePriorityQueue>(sortedStates);
             }
             
+            std::shared_ptr<StatePriorityQueue> createStatePriorityQueue(std::vector<storm::storage::sparse::state_type> const& states) {
+                return std::make_shared<StaticStatePriorityQueue>(states);
+            }
+            
             template uint_fast64_t estimateComplexity(double const& value);
             template std::shared_ptr<StatePriorityQueue> createStatePriorityQueue(boost::optional<std::vector<uint_fast64_t>> const& distanceBasedStatePriorities, storm::storage::FlexibleSparseMatrix<double> const& transitionMatrix, storm::storage::FlexibleSparseMatrix<double> const& backwardTransitions, std::vector<double>& oneStepProbabilities, storm::storage::BitVector const& states);
             template uint_fast64_t computeStatePenalty(storm::storage::sparse::state_type const& state, storm::storage::FlexibleSparseMatrix<double> const& transitionMatrix, storm::storage::FlexibleSparseMatrix<double> const& backwardTransitions, std::vector<double> const& oneStepProbabilities);
             template uint_fast64_t computeStatePenaltyRegularExpression(storm::storage::sparse::state_type const& state, storm::storage::FlexibleSparseMatrix<double> const& transitionMatrix, storm::storage::FlexibleSparseMatrix<double> const& backwardTransitions, std::vector<double> const& oneStepProbabilities);
 
+            template uint_fast64_t estimateComplexity(storm::RationalNumber const& value);
             
 #ifdef STORM_HAVE_CARL
             template uint_fast64_t estimateComplexity(storm::RationalFunction const& value);

src/utility/stateelimination.h

@@ -53,7 +53,8 @@ namespace storm {
             template<typename ValueType>
             std::shared_ptr<StatePriorityQueue> createStatePriorityQueue(boost::optional<std::vector<uint_fast64_t>> const& stateDistances, storm::storage::FlexibleSparseMatrix<ValueType> const& transitionMatrix, storm::storage::FlexibleSparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector const& states);
             
-            std::shared_ptr<StatePriorityQueue> createNaivePriorityQueue(storm::storage::BitVector const& states);
+            std::shared_ptr<StatePriorityQueue> createStatePriorityQueue(storm::storage::BitVector const& states);
+            std::shared_ptr<StatePriorityQueue> createStatePriorityQueue(std::vector<storm::storage::sparse::state_type> const& states);
             
         }
     }

test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.cpp

@@ -0,0 +1,321 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/parser/FormulaParser.h"
+#include "src/logic/Formulas.h"
+#include "src/utility/solver.h"
+#include "src/models/sparse/StandardRewardModel.h"
+#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
+#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
+#include "src/settings/SettingsManager.h"
+#include "src/settings/modules/GeneralSettings.h"
+#include "src/settings/modules/EigenEquationSolverSettings.h"
+#include "src/settings/modules/NativeEquationSolverSettings.h"
+#include "src/settings/SettingMemento.h"
+#include "src/parser/AutoParser.h"
+#include "src/parser/PrismParser.h"
+#include "src/builder/ExplicitPrismModelBuilder.h"
+
+TEST(EigenDtmcPrctlModelCheckerTest, Die) {
+    std::shared_ptr<storm::models::sparse::Model<double>> abstractModel = storm::parser::AutoParser<>::parseModel(STORM_CPP_BASE_PATH "/examples/dtmc/die/die.tra", STORM_CPP_BASE_PATH "/examples/dtmc/die/die.lab", "", STORM_CPP_BASE_PATH "/examples/dtmc/die/die.coin_flips.trans.rew");
+    
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser formulaParser;
+    
+    ASSERT_EQ(abstractModel->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc = abstractModel->as<storm::models::sparse::Dtmc<double>>();
+    
+    ASSERT_EQ(dtmc->getNumberOfStates(), 13ull);
+    ASSERT_EQ(dtmc->getNumberOfTransitions(), 20ull);
+    
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+    
+    std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(1.0 / 6.0, quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    
+    formula = formulaParser.parseSingleFormulaFromString("P=? [F \"two\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(1.0 / 6.0, quantitativeResult2[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    
+    formula = formulaParser.parseSingleFormulaFromString("P=? [F \"three\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(1.0 / 6.0, quantitativeResult3[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    
+    formula = formulaParser.parseSingleFormulaFromString("R=? [F \"done\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult4 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(11.0 / 3.0, quantitativeResult4[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+}
+
+TEST(EigenDtmcPrctlModelCheckerTest, Crowds) {
+    std::shared_ptr<storm::models::sparse::Model<double>> abstractModel = storm::parser::AutoParser<>::parseModel(STORM_CPP_BASE_PATH "/examples/dtmc/crowds/crowds5_5.tra", STORM_CPP_BASE_PATH "/examples/dtmc/crowds/crowds5_5.lab", "", "");
+    
+    ASSERT_EQ(abstractModel->getType(), storm::models::ModelType::Dtmc);
+    
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser formulaParser;
+    
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc = abstractModel->as<storm::models::sparse::Dtmc<double>>();
+    
+    ASSERT_EQ(8607ull, dtmc->getNumberOfStates());
+    ASSERT_EQ(15113ull, dtmc->getNumberOfTransitions());
+    
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+    
+    std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observe0Greater1\"]");
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(0.3328800375801578281, quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    
+    formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observeIGreater1\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(0.1522194965, quantitativeResult2[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    
+    formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observeOnlyTrueSender\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(0.32153724292835045, quantitativeResult3[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+}
+
+TEST(EigenDtmcPrctlModelCheckerTest, SynchronousLeader) {
+    std::shared_ptr<storm::models::sparse::Model<double>> abstractModel = storm::parser::AutoParser<>::parseModel(STORM_CPP_BASE_PATH "/examples/dtmc/synchronous_leader/leader4_8.tra", STORM_CPP_BASE_PATH "/examples/dtmc/synchronous_leader/leader4_8.lab", "", STORM_CPP_BASE_PATH "/examples/dtmc/synchronous_leader/leader4_8.pick.trans.rew");
+    
+    ASSERT_EQ(abstractModel->getType(), storm::models::ModelType::Dtmc);
+    
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser formulaParser;
+    
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc = abstractModel->as<storm::models::sparse::Dtmc<double>>();
+    
+    ASSERT_EQ(12400ull, dtmc->getNumberOfStates());
+    ASSERT_EQ(16495ull, dtmc->getNumberOfTransitions());
+    
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+    
+    std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"elected\"]");
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(1.0, quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    
+    formula = formulaParser.parseSingleFormulaFromString("P=? [F<=20 \"elected\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(0.9999965911265462636, quantitativeResult2[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    
+    formula = formulaParser.parseSingleFormulaFromString("R=? [F \"elected\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
+    
+    EXPECT_NEAR(1.0448979591836789, quantitativeResult3[0], storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+}
+
+TEST(EigenDtmcPrctlModelCheckerTest, LRASingleBscc) {
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc;
+    
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser formulaParser;
+    
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 2);
+        matrixBuilder.addNextValue(0, 1, 1.);
+        matrixBuilder.addNextValue(1, 0, 1.);
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(2);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 1);
+        
+        dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
+        
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        
+        std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+    }
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 4);
+        matrixBuilder.addNextValue(0, 0, .5);
+        matrixBuilder.addNextValue(0, 1, .5);
+        matrixBuilder.addNextValue(1, 0, .5);
+        matrixBuilder.addNextValue(1, 1, .5);
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(2);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 1);
+        
+        dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
+        
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        
+        std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+    }
+    
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(3, 3, 3);
+        matrixBuilder.addNextValue(0, 1, 1);
+        matrixBuilder.addNextValue(1, 2, 1);
+        matrixBuilder.addNextValue(2, 0, 1);
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(3);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 2);
+        
+        dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
+        
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+        
+        std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+    }
+}
+
+TEST(EigenDtmcPrctlModelCheckerTest, LRA) {
+    storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc;
+    
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser formulaParser;
+    
+    {
+        matrixBuilder = storm::storage::SparseMatrixBuilder<double>(15, 15, 20, true);
+        matrixBuilder.addNextValue(0, 1, 1);
+        matrixBuilder.addNextValue(1, 4, 0.7);
+        matrixBuilder.addNextValue(1, 6, 0.3);
+        matrixBuilder.addNextValue(2, 0, 1);
+        
+        matrixBuilder.addNextValue(3, 5, 0.8);
+        matrixBuilder.addNextValue(3, 9, 0.2);
+        matrixBuilder.addNextValue(4, 3, 1);
+        matrixBuilder.addNextValue(5, 3, 1);
+        
+        matrixBuilder.addNextValue(6, 7, 1);
+        matrixBuilder.addNextValue(7, 8, 1);
+        matrixBuilder.addNextValue(8, 6, 1);
+        
+        matrixBuilder.addNextValue(9, 10, 1);
+        matrixBuilder.addNextValue(10, 9, 1);
+        matrixBuilder.addNextValue(11, 9, 1);
+        
+        matrixBuilder.addNextValue(12, 5, 0.4);
+        matrixBuilder.addNextValue(12, 8, 0.3);
+        matrixBuilder.addNextValue(12, 11, 0.3);
+        
+        matrixBuilder.addNextValue(13, 7, 0.7);
+        matrixBuilder.addNextValue(13, 12, 0.3);
+        
+        matrixBuilder.addNextValue(14, 12, 1);
+        
+        storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+        
+        storm::models::sparse::StateLabeling ap(15);
+        ap.addLabel("a");
+        ap.addLabelToState("a", 1);
+        ap.addLabelToState("a", 4);
+        ap.addLabelToState("a", 5);
+        ap.addLabelToState("a", 7);
+        ap.addLabelToState("a", 11);
+        ap.addLabelToState("a", 13);
+        ap.addLabelToState("a", 14);
+        
+        dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
+        
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+        
+        std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
+        
+        std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+        storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+        
+        EXPECT_NEAR(0.3 / 3., quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(0.0, quantitativeResult1[3], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(1. / 3., quantitativeResult1[6], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(0.0, quantitativeResult1[9], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(0.3 / 3., quantitativeResult1[12], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(.79 / 3., quantitativeResult1[13], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+        EXPECT_NEAR(0.3 / 3., quantitativeResult1[14], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+    }
+}
+
+TEST(EigenDtmcPrctlModelCheckerTest, Conditional) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/modelchecker/test_conditional.pm");
+    
+    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::builder::ExplicitPrismModelBuilder<double>(program).translate();
+    ASSERT_TRUE(model->getType() == storm::models::ModelType::Dtmc);
+    ASSERT_EQ(4ul, model->getNumberOfStates());
+    ASSERT_EQ(5ul, model->getNumberOfTransitions());
+    
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc = model->as<storm::models::sparse::Dtmc<double>>();
+    
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+    
+    // A parser that we use for conveniently constructing the formulas.
+    storm::parser::FormulaParser formulaParser;
+    std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"target\"]");
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(0.5, quantitativeResult1[0], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+    
+    formula = formulaParser.parseSingleFormulaFromString("P=? [F \"target\" || F \"condition\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(storm::utility::one<double>(), quantitativeResult2[0], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+    
+    formula = formulaParser.parseSingleFormulaFromString("R=? [F \"target\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_EQ(storm::utility::infinity<double>(), quantitativeResult3[0]);
+    
+    formula = formulaParser.parseSingleFormulaFromString("R=? [F \"target\" || F \"condition\"]");
+    
+    result = checker.check(*formula);
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult4 = result->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(storm::utility::one<double>(), quantitativeResult4[0], storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
+}

test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.hpp

@@ -0,0 +1,14 @@
+//
+//  EigenDtmcPrctlModelCheckerTest.hpp
+//  storm
+//
+//  Created by Christian Dehnert on 19.06.16.
+//
+//
+
+#ifndef EigenDtmcPrctlModelCheckerTest_hpp
+#define EigenDtmcPrctlModelCheckerTest_hpp
+
+#include <stdio.h>
+
+#endif /* EigenDtmcPrctlModelCheckerTest_hpp */

test/functional/solver/EigenLinearEquationSolverTest.cpp

@@ -4,7 +4,8 @@
 #include "src/solver/EigenLinearEquationSolver.h"
 #include "src/settings/SettingsManager.h"
 
-#include "src/settings/modules/GmmxxEquationSolverSettings.h"
+#include "src/utility/constants.h"
+#include "src/settings/modules/EigenEquationSolverSettings.h"
 
 TEST(EigenLinearEquationSolver, SolveWithStandardOptions) {
     ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
@@ -29,9 +30,9 @@ TEST(EigenLinearEquationSolver, SolveWithStandardOptions) {
     
     storm::solver::EigenLinearEquationSolver<double> solver(A);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
-    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
 }
 
 TEST(EigenLinearEquationSolver, SparseLU) {
@@ -57,9 +58,65 @@ TEST(EigenLinearEquationSolver, SparseLU) {
     
     storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::SparseLU, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
-    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+}
+
+TEST(EigenLinearEquationSolver, SparseLU_RationalNumber) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<storm::RationalNumber> builder);
+    storm::storage::SparseMatrixBuilder<storm::RationalNumber> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, 2));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 4));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, -2));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, 4));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, 3));
+    
+    storm::storage::SparseMatrix<storm::RationalNumber> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<storm::RationalNumber> x(3);
+    std::vector<storm::RationalNumber> b = {16, -4, -7};
+    
+    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<storm::RationalNumber> solver(A, storm::solver::EigenLinearEquationSolver<storm::RationalNumber>::SolutionMethod::SparseLU));
+    
+    storm::solver::EigenLinearEquationSolver<storm::RationalNumber> solver(A, storm::solver::EigenLinearEquationSolver<storm::RationalNumber>::SolutionMethod::SparseLU);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
+    ASSERT_TRUE(storm::utility::isOne(x[0]));
+    ASSERT_TRUE(x[1] == 3);
+    ASSERT_TRUE(x[2] == -1);
+}
+
+TEST(EigenLinearEquationSolver, SparseLU_RationalFunction) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<storm::RationalFunction> builder);
+    storm::storage::SparseMatrixBuilder<storm::RationalFunction> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, storm::RationalFunction(2)));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, storm::RationalFunction(4)));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, storm::RationalFunction(-2)));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, storm::RationalFunction(4)));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, storm::RationalFunction(-1)));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, storm::RationalFunction(5)));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, storm::RationalFunction(-1)));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, storm::RationalFunction(-1)));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, storm::RationalFunction(3)));
+    
+    storm::storage::SparseMatrix<storm::RationalFunction> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<storm::RationalFunction> x(3);
+    std::vector<storm::RationalFunction> b = {storm::RationalFunction(16), storm::RationalFunction(-4), storm::RationalFunction(-7)};
+    
+    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<storm::RationalFunction> solver(A, storm::solver::EigenLinearEquationSolver<storm::RationalFunction>::SolutionMethod::SparseLU));
+    
+    storm::solver::EigenLinearEquationSolver<storm::RationalFunction> solver(A, storm::solver::EigenLinearEquationSolver<storm::RationalFunction>::SolutionMethod::SparseLU);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
+    ASSERT_TRUE(storm::utility::isOne(x[0]));
+    ASSERT_TRUE(x[1] == storm::RationalFunction(3));
+    ASSERT_TRUE(x[2] == storm::RationalFunction(-1));
 }
 
 TEST(EigenLinearEquationSolver, BiCGSTAB) {
@@ -85,9 +142,9 @@ TEST(EigenLinearEquationSolver, BiCGSTAB) {
     
     storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
-    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
 }
 
 TEST(EigenLinearEquationSolver, BiCGSTAB_Ilu) {
@@ -113,9 +170,9 @@ TEST(EigenLinearEquationSolver, BiCGSTAB_Ilu) {
     
     storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Ilu);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
-    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
 }
 
 TEST(EigenLinearEquationSolver, BiCGSTAB_Diagonal) {
@@ -141,9 +198,9 @@ TEST(EigenLinearEquationSolver, BiCGSTAB_Diagonal) {
     
     storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Diagonal);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
-    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
-    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
 }
 
 TEST(EigenLinearEquationSolver, MatrixVectorMultplication) {
@@ -166,5 +223,5 @@ TEST(EigenLinearEquationSolver, MatrixVectorMultplication) {
     
     storm::solver::EigenLinearEquationSolver<double> solver(A);
     ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(x, nullptr, 4));
-    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
 }

test/functional/solver/EliminationLinearEquationSolverTest.cpp

@@ -0,0 +1,58 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/solver/EliminationLinearEquationSolver.h"
+#include "src/settings/SettingsManager.h"
+
+#include "src/settings/modules/GmmxxEquationSolverSettings.h"
+
+TEST(EliminationLinearEquationSolver, Solve) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, 2));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 4));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, -2));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, 4));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, 3));
+    
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(3);
+    std::vector<double> b = {16, -4, -7};
+    
+    ASSERT_NO_THROW(storm::solver::EliminationLinearEquationSolver<double> solver(A));
+    
+    storm::solver::EliminationLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
+    ASSERT_LT(std::abs(x[0] - 1), 1e-15);
+    ASSERT_LT(std::abs(x[1] - 3), 1e-15);
+    ASSERT_LT(std::abs(x[2] - (-1)), 1e-15);
+}
+
+TEST(EliminationLinearEquationSolver, MatrixVectorMultplication) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(0, 4, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(1, 4, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 3, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 4, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(3, 4, 1));
+    ASSERT_NO_THROW(builder.addNextValue(4, 4, 1));
+    
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(5);
+    x[4] = 1;
+    
+    storm::solver::EliminationLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(x, nullptr, 4));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+}