Refactored linear equation solvers and nondeterministic linear equation solvers. Added functional tests for both.

Former-commit-id: 0abb11828a

resources/3rdparty/gmm-4.2/include/gmm/gmm_blas.h

@@ -38,7 +38,7 @@
 #ifndef GMM_BLAS_H__
 #define GMM_BLAS_H__
 
-// This Version of GMM was modified for StoRM
+// This Version of GMM was modified for StoRM.
 // To detect whether the usage of TBB is possible, this include is neccessary
 #include "storm-config.h"
 
@@ -404,64 +404,6 @@ namespace gmm {
     for (; it != ite; ++it, ++it2) res += (*it) * (*it2);
     return res;
   }
-  
-#ifdef STORM_HAVE_INTELTBB
-	/* Official Intel Hint on blocked_range vs. linear iterators: http://software.intel.com/en-us/forums/topic/289505
-
-	 */
-	template <typename IT1>
-class forward_range {
-    IT1 my_begin;
-    IT1 my_end;
-    size_t my_size;
-public:
-    IT1 begin() const {return my_begin;}
-    IT1 end() const {return my_end;}
-    bool empty() const {return my_begin==my_end;}
-    bool is_divisible() const {return my_size>1;}
-    forward_range( IT1 first, IT1 last, size_t size ) : my_begin(first), my_end(last), my_size(size) {
-        assert( size==size_t(std::distance( first,last )));
-    }
-	forward_range( IT1 first, IT1 last) : my_begin(first), my_end(last) {
-		my_size = std::distance( first,last );
-    }
-    forward_range( forward_range& r, tbb::split ) {
-        size_t h = r.my_size/2;
-        my_end = r.my_end;
-        my_begin = r.my_begin;
-        std::advance( my_begin, h ); // Might be scaling issue
-        my_size = r.my_size-h;
-        r.my_end = my_begin;
-        r.my_size = h;
-    }
-};
-
-	template <typename IT1, typename V>
-	class tbbHelper_vect_sp_sparse {
-		V const* my_v;
-		public:
-			 typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
-				typename linalg_traits<V>::value_type>::T my_sum; 
-			void operator()( const forward_range<IT1>& r ) {
-				V const* v = my_v;
-				typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
-				typename linalg_traits<V>::value_type>::T sum = my_sum;
-				IT1 end = r.end();
-				for( IT1 i=r.begin(); i!=end; ++i) { 
-					sum += (*i) * v->at(i.index());
-				}
-				my_sum = sum;    
-			}
- 
-			tbbHelper_vect_sp_sparse( tbbHelper_vect_sp_sparse& x, tbb::split ) : my_v(x.my_v), my_sum(0) {}
- 
-			void join( const tbbHelper_vect_sp_sparse& y ) {my_sum+=y.my_sum;}
-             
-			tbbHelper_vect_sp_sparse(V const* v) :
-				my_v(v), my_sum(0)
-			{}
-	};
-#endif
 
   template <typename IT1, typename V> inline
     typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
@@ -469,14 +411,7 @@ public:
     vect_sp_sparse_(IT1 it, IT1 ite, const V &v) {
       typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type,
 	typename linalg_traits<V>::value_type>::T res(0);
-#if defined(STORM_HAVE_INTELTBB) && defined(STORM_USE_TBB_FOR_INNER)
-	  // This is almost never an efficent way, only if there are _many_ states
-	  tbbHelper_vect_sp_sparse<IT1, V> tbbHelper(&v);
-	  tbb::parallel_reduce(forward_range<IT1>(it, ite), tbbHelper);
-	  res = tbbHelper.my_sum;
-#else
     for (; it != ite; ++it) res += (*it) * v[it.index()];
-#endif
     return res;
   }
 

src/adapters/GmmxxAdapter.h

@@ -67,7 +67,6 @@ public:
 	template<class T>
 	static std::unique_ptr<gmm::csr_matrix<T>> toGmmxxSparseMatrix(storm::storage::SparseMatrix<T>&& matrix) {
 		uint_fast64_t realNonZeros = matrix.getEntryCount();
-        std::cout << "here?!" << std::endl;
 		LOG4CPLUS_DEBUG(logger, "Converting matrix with " << realNonZeros << " non-zeros to gmm++ format.");
 
 		// Prepare the resulting matrix.
@@ -77,11 +76,7 @@ public:
         typedef std::vector<IND_TYPE> vectorType_ull_t;
         typedef std::vector<T> vectorType_T_t; 
 
-        // Move Row Indications
-        result->jc.~vectorType_ull_t(); // Call Destructor inplace
-		new (&result->jc) vectorType_ull_t(std::move(*storm::utility::ConversionHelper::toUnsignedLongLong(&matrix.rowIndications)));
-        
-        // Copy columns and values.
+        // Copy columns and values. It is absolutely necessary to do so before moving the row indications vector.
         std::vector<T> values;
         values.reserve(matrix.getEntryCount());
         std::vector<uint_fast64_t> columns;
@@ -91,6 +86,9 @@ public:
             columns.emplace_back(entry.first);
             values.emplace_back(entry.second);
         }
+
+        // Move Row Indications
+        result->jc = std::move(matrix.rowIndications);
         
         std::swap(result->ir, columns);
         std::swap(result->pr, values);

src/counterexamples/PathBasedSubsystemGenerator.h

@@ -483,8 +483,6 @@ public:
 
 		//At last compensate for the distance between init and source state
 		probability = itSearch ? probability : probability / itDistances[bestIndex].second;
-
-		LOG4CPLUS_DEBUG(logger, "Found path: " << shortestPath);
 	}
 
 private:

src/solver/AbstractNondeterministicLinearEquationSolver.h

@@ -1,148 +1,64 @@
 #ifndef STORM_SOLVER_ABSTRACTNONDETERMINISTICLINEAREQUATIONSOLVER_H_
 #define STORM_SOLVER_ABSTRACTNONDETERMINISTICLINEAREQUATIONSOLVER_H_
 
-#include "src/storage/SparseMatrix.h"
-#include "src/utility/vector.h"
-#include "src/settings/Settings.h"
-
 #include <vector>
 
+#include "src/storage/SparseMatrix.h"
+
 namespace storm {
     namespace solver {
         
-        template<class Type>
+        /*!
+         * A class that represents an abstract nondeterministic linear equation solver. In addition to solving linear
+         * equation systems involving min/max operators, .
+         */
+        template<class ValueType>
         class AbstractNondeterministicLinearEquationSolver {
         public:
-            AbstractNondeterministicLinearEquationSolver() {
-                storm::settings::Settings* s = storm::settings::Settings::getInstance();
-                precision = s->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
-                maxIterations = s->getOptionByLongName("maxiter").getArgument(0).getValueAsUnsignedInteger();
-                relative = !s->isSet("absolute");
-            }
-            
-            AbstractNondeterministicLinearEquationSolver(double precision, uint_fast64_t maxIterations, bool relative) : precision(precision), maxIterations(maxIterations), relative(relative) {
-                // Intentionally left empty.
-            }
-            
-            virtual AbstractNondeterministicLinearEquationSolver<Type>* clone() const {
-                return new AbstractNondeterministicLinearEquationSolver<Type>(this->precision, this->maxIterations, this->relative);
-            }
-            
             /*!
-             * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes x[i+1] = A*x[i] + b
-             * until x[n], where x[0] = x.
+             * Makes a copy of the nondeterministic linear equation solver.
              *
-             * @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise.
-             * @param A The matrix that is to be multiplied against the vector.
-             * @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter,
-             * i.e. after the method returns, this vector will contain the computed values.
-             * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix.
-             * @param b If not null, this vector is being added to the result after each matrix-vector multiplication.
-             * @param n Specifies the number of iterations the matrix-vector multiplication is performed.
-             * @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector.
+             * @return A pointer to a copy of the nondeterministic linear equation solver.
              */
-            virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const {
-                // Create vector for result of multiplication, which is reduced to the result vector after
-                // each multiplication.
-                std::vector<Type> multiplyResult(A.getRowCount());
-                
-                // Now perform matrix-vector multiplication as long as we meet the bound of the formula.
-                for (uint_fast64_t i = 0; i < n; ++i) {
-                    A.multiplyWithVector(x, multiplyResult);
-                    
-                    // Add b if it is non-null.
-                    if (b != nullptr) {
-                        storm::utility::vector::addVectorsInPlace(multiplyResult, *b);
-                    }
-                    
-                    // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
-                    // element of the min/max operator stack.
-                    if (minimize) {
-                        storm::utility::vector::reduceVectorMin(multiplyResult, x, nondeterministicChoiceIndices);
-                    } else {
-                        storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices);
-                    }
-                }
-            }
+            virtual AbstractNondeterministicLinearEquationSolver<ValueType>* clone() const = 0;
             
             /*!
-             * Solves the equation system A*x = b given by the parameters.
+             * Solves the equation system x = min/max(A*x + b) given by the parameters.
              *
-             * @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise.
+             * @param minimize If set, all the value of a group of rows is the taken as the minimum over all rows and as
+             * the maximum otherwise.
              * @param A The matrix specifying the coefficients of the linear equations.
-             * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but
-             * may be ignored.
-             * @param b The right-hand side of the equation system.
-             * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix.
-             * @param takenChoices If not null, this vector will be filled with the nondeterministic choices taken by the states
-             * to achieve the solution of the equation system. This assumes that the given vector has at least as many elements
-             * as there are states in the MDP.
-             * @returns The solution vector x of the system of linear equations as the content of the parameter x.
+             * @param x The solution vector x. The initial values of x represent a guess of the real values to the
+             * solver, but may be ignored.
+             * @param b The vector to add after matrix-vector multiplication.
+             * @param rowGroupIndices A vector indicating which rows of the matrix belong to one group.
+             * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
+             * vector must be equal to the number of rows of A.
+             * @param newX If non-null, this memory is used as a scratch memory. If given, the length of this
+             * vector must be equal to the length of the vector x (and thus to the number of columns of A).
+             * @return The solution vector x of the system of linear equations as the content of the parameter x.
              */
-            virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* multiplyResult = nullptr, std::vector<Type>* newX = nullptr) const {
-                
-                // Set up the environment for the power method.
-                bool multiplyResultMemoryProvided = true;
-                if (multiplyResult == nullptr) {
-                    multiplyResult = new std::vector<Type>(A.getRowCount());
-                    multiplyResultMemoryProvided = false;
-                }
-                std::vector<Type>* currentX = &x;
-                bool xMemoryProvided = true;
-                if (newX == nullptr) {
-                    newX = new std::vector<Type>(x.size());
-                    xMemoryProvided = false;
-                }
-                std::vector<Type>* swap = nullptr;
-                uint_fast64_t iterations = 0;
-                bool converged = false;
-
-                // Proceed with the iterations as long as the method did not converge or reach the
-                // user-specified maximum number of iterations.
-                while (!converged && iterations < maxIterations) {
-                    // Compute x' = A*x + b.
-                    A.multiplyWithVector(*currentX, *multiplyResult);
-                    storm::utility::vector::addVectorsInPlace(*multiplyResult, b);
-                    
-                    // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
-                    // element of the min/max operator stack.
-                    if (minimize) {
-                        storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, nondeterministicChoiceIndices);
-                    } else {
-                        storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, nondeterministicChoiceIndices);
-                    }
-                    
-                    // Determine whether the method converged.
-                    converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative);
-                    
-                    // Update environment variables.
-                    swap = currentX;
-                    currentX = newX;
-                    newX = swap;
-                    ++iterations;
-                }
-
-                // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
-                // is currently stored in currentX, but x is the output vector.
-                if (iterations % 2 == 1) {
-                    std::swap(x, *currentX);
-                    if (!xMemoryProvided) {
-                        delete currentX;
-                    }
-                } else if (!xMemoryProvided) {
-                    delete newX;
-                }
-                
-                if (!multiplyResultMemoryProvided) {
-                    delete multiplyResult;
-                }
-            }
+            virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<uint_fast64_t> const& rowGroupIndices, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const = 0;
             
-        protected:
-            double precision;
-            uint_fast64_t maxIterations;
-            bool relative;
-
+            /*!
+             * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes
+             * x[i+1] = min/max(A*x[i] + b) until x[n], where x[0] = x. After each multiplication and addition, the
+             * minimal/maximal value out of each row group is selected to reduce the resulting vector to obtain the
+             * vector for the next iteration.
+             *
+             * @param minimize If set, all the value of a group of rows is the taken as the minimum over all rows and as
+             * the maximum otherwise.
+             * @param A The matrix that is to be multiplied with the vector.
+             * @param x The initial vector that is to be multiplied with the matrix. This is also the output parameter,
+             * i.e. after the method returns, this vector will contain the computed values.
+             * @param rowGroupIndices A vector indicating which rows of the matrix belong to one group.
+             * @param b If not null, this vector is added after each multiplication.
+             * @param n Specifies the number of iterations the matrix-vector multiplication is performed.
+             * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
+             * vector must be equal to the number of rows of A.
+             * @return The result of the repeated matrix-vector multiplication as the content of the vector x.
+             */
+            virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<uint_fast64_t> const& rowGroupIndices, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const = 0;
         };
         
     } // namespace solver

src/solver/GmmxxLinearEquationSolver.cpp

@@ -45,7 +45,7 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(SolutionMethod method, double precision, bool relative, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, uint_fast64_t restart) : method(method), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner), restart(restart) {
+        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, bool relative, uint_fast64_t restart) : method(method), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner), restart(restart) {
             // Intentionally left empty.
         }
         
@@ -150,12 +150,12 @@ namespace storm {
         void GmmxxLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
             // Transform the transition probability A to the gmm++ format to use its arithmetic.
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
-            
+
             // Set up some temporary variables so that we can just swap pointers instead of copying the result after
             // each iteration.
             std::vector<ValueType>* swap = nullptr;
             std::vector<ValueType>* currentX = &x;
-            
+
             bool multiplyResultProvided = true;
             std::vector<ValueType>* nextX = multiplyResult;
             if (nextX == nullptr) {
@@ -182,7 +182,7 @@ namespace storm {
             }
             
             // If the vector for the temporary multiplication result was not provided, we need to delete it.
-            if (multiplyResultProvided) {
+            if (!multiplyResultProvided) {
                 delete copyX;
             }
         }
@@ -196,7 +196,7 @@ namespace storm {
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmDinv = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.second));
             // Convert the LU matrix to gmm++'s format.
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.first));
-            
+        
             // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
             bool multiplyResultProvided = true;
             std::vector<ValueType>* nextX = multiplyResult;
@@ -220,12 +220,12 @@ namespace storm {
                 gmm::add(b, gmm::scaled(tmpX, -storm::utility::constantOne<ValueType>()), tmpX);
                 gmm::mult(*gmmDinv, tmpX, *nextX);
                 
-                // Swap the two pointers as a preparation for the next iteration.
-                std::swap(nextX, currentX);
-                
                 // Now check if the process already converged within our precision.
                 converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, precision, relative);
-                
+
+                // Swap the two pointers as a preparation for the next iteration.
+                std::swap(nextX, currentX);
+
                 // Increase iteration count so we can abort if convergence is too slow.
                 ++iterationCount;
             }

src/solver/GmmxxLinearEquationSolver.h

@@ -27,14 +27,14 @@ namespace storm {
              *
              * @param method The method to use for linear equation solving.
              * @param precision The precision to use for convergence detection.
-             * @param relative If set, the relative error rather than the absolute error is considered for convergence
-             * detection.
              * @param maximalNumberOfIterations The maximal number of iterations do perform before iteration is aborted.
              * @param preconditioner The preconditioner to use when solving linear equation systems.
+             * @param relative If set, the relative error rather than the absolute error is considered for convergence
+             * detection.
              * @param restart An optional argument that specifies after how many iterations restarted methods are
              * supposed to actually to a restart.
              */
-            GmmxxLinearEquationSolver(SolutionMethod method, double precision, bool relative, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, uint_fast64_t restart = 0);
+            GmmxxLinearEquationSolver(SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, bool relative = true, uint_fast64_t restart = 0);
             
             /*!
              * Constructs a linear equation solver with parameters being set according to the settings object.

src/solver/GmmxxNondeterministicLinearEquationSolver.cpp

@@ -1,11 +1,144 @@
+#include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
 #include "src/settings/Settings.h"
+#include "src/adapters/GmmxxAdapter.h"
+#include "src/utility/vector.h"
 
 bool GmmxxNondeterministicLinearEquationSolverOptionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* instance) -> bool {
 	instance->addOption(storm::settings::OptionBuilder("GmmxxNondeterminsticLinearEquationSolver", "maxiter", "i", "The maximal number of iterations to perform before iterative solving is aborted.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(10000).build()).build());
     
-	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "precision", "", "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-6).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+	instance->addOption(storm::settings::OptionBuilder("GmmxxNondeterminsticLinearEquationSolver", "precision", "", "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-6).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
     
-	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "absolute", "", "Whether the relative or the absolute error is considered for deciding convergence.").build());
+	instance->addOption(storm::settings::OptionBuilder("GmmxxNondeterminsticLinearEquationSolver", "absolute", "", "Whether the relative or the absolute error is considered for deciding convergence.").build());
     
 	return true;
-});
+});
+
+namespace storm {
+    namespace solver {
+        
+        template<typename ValueType>
+        GmmxxNondeterministicLinearEquationSolver<ValueType>::GmmxxNondeterministicLinearEquationSolver() {
+            // Get the settings object to customize solving.
+            storm::settings::Settings* settings = storm::settings::Settings::getInstance();
+            
+            // Get appropriate settings.
+            maximalNumberOfIterations = settings->getOptionByLongName("maxiter").getArgument(0).getValueAsUnsignedInteger();
+            precision = settings->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
+            relative = settings->isSet("absolute");
+        }
+        
+        template<typename ValueType>
+        GmmxxNondeterministicLinearEquationSolver<ValueType>::GmmxxNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
+            // Intentionally left empty.
+        }
+
+        
+        template<typename ValueType>
+        AbstractNondeterministicLinearEquationSolver<ValueType>* GmmxxNondeterministicLinearEquationSolver<ValueType>::clone() const {
+            return new GmmxxNondeterministicLinearEquationSolver<ValueType>(*this);
+        }
+        
+        template<typename ValueType>
+        void GmmxxNondeterministicLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+            // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
+            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
+            
+            // Set up the environment for the power method. If scratch memory was not provided, we need to create it.
+            bool multiplyResultMemoryProvided = true;
+            if (multiplyResult == nullptr) {
+                multiplyResult = new std::vector<ValueType>(A.getRowCount());
+                multiplyResultMemoryProvided = false;
+            }
+            
+            std::vector<ValueType>* currentX = &x;
+            bool xMemoryProvided = true;
+            if (newX == nullptr) {
+                newX = new std::vector<ValueType>(x.size());
+                xMemoryProvided = false;
+            }
+            std::vector<ValueType>* swap = nullptr;
+            uint_fast64_t iterations = 0;
+            bool converged = false;
+            
+            // Keep track of which of the vectors for x is the auxiliary copy.
+            std::vector<ValueType>* copyX = newX;
+
+            // Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number
+            // of iterations.
+            while (!converged && iterations < maximalNumberOfIterations) {
+                // Compute x' = A*x + b.
+                gmm::mult(*gmmxxMatrix, *currentX, *multiplyResult);
+                gmm::add(b, *multiplyResult);
+                
+                // Reduce the vector x by applying min/max over all nondeterministic choices.
+                if (minimize) {
+                    storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, nondeterministicChoiceIndices);
+                } else {
+                    storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, nondeterministicChoiceIndices);
+                }
+                
+                // Determine whether the method converged.
+                converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->precision, this->relative);
+                
+                // Update environment variables.
+                std::swap(currentX, newX);
+                ++iterations;
+            }
+            
+            // Check if the solver converged and issue a warning otherwise.
+            if (converged) {
+                LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
+            } else {
+                LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations.");
+            }
+            
+            // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
+            // is currently stored in currentX, but x is the output vector.
+            if (currentX == copyX) {
+                std::swap(x, *currentX);
+            }
+            
+            if (!xMemoryProvided) {
+                delete copyX;
+            }
+            
+            if (!multiplyResultMemoryProvided) {
+                delete multiplyResult;
+            }
+        }
+        
+        template<typename ValueType>
+        void GmmxxNondeterministicLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+            // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
+            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
+            
+            bool multiplyResultMemoryProvided = true;
+            if (multiplyResult == nullptr) {
+                multiplyResult = new std::vector<ValueType>(A.getRowCount());
+                multiplyResultMemoryProvided = false;
+            }
+            
+            // Now perform matrix-vector multiplication as long as we meet the bound of the formula.
+            for (uint_fast64_t i = 0; i < n; ++i) {
+                gmm::mult(*gmmxxMatrix, x, *multiplyResult);
+                
+                if (b != nullptr) {
+                    gmm::add(*b, *multiplyResult);
+                }
+                
+                if (minimize) {
+                    storm::utility::vector::reduceVectorMin(*multiplyResult, x, nondeterministicChoiceIndices);
+                } else {
+                    storm::utility::vector::reduceVectorMax(*multiplyResult, x, nondeterministicChoiceIndices);
+                }
+            }
+            
+            if (!multiplyResultMemoryProvided) {
+                delete multiplyResult;
+            }
+        }
+
+        // Explicitly instantiate the solver.
+        template class GmmxxNondeterministicLinearEquationSolver<double>;
+    } // namespace solver
+} // namespace storm

src/solver/GmmxxNondeterministicLinearEquationSolver.h

@@ -1,13 +1,7 @@
 #ifndef STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_
 #define STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_
 
-#include "AbstractLinearEquationSolver.h"
-#include "src/adapters/GmmxxAdapter.h"
-#include "src/storage/SparseMatrix.h"
-
-#include "gmm/gmm_matrix.h"
-
-#include <vector>
+#include "src/solver/AbstractNondeterministicLinearEquationSolver.h"
 
 namespace storm {
     namespace solver {
@@ -15,105 +9,38 @@ namespace storm {
         template<class Type>
         class GmmxxNondeterministicLinearEquationSolver : public AbstractNondeterministicLinearEquationSolver<Type> {
         public:
+            /*!
+             * Constructs a nondeterministic linear equation solver with parameters being set according to the settings
+             * object.
+             */
+            GmmxxNondeterministicLinearEquationSolver();
             
-            virtual AbstractNondeterministicLinearEquationSolver<Type>* clone() const {
-                return new GmmxxNondeterministicLinearEquationSolver<Type>();
-            }
+            /*!
+             * Constructs a nondeterminstic linear equation solver with the given parameters.
+             *
+             * @param precision The precision to use for convergence detection.
+             * @param maximalNumberOfIterations The maximal number of iterations do perform before iteration is aborted.
+             * @param relative If set, the relative error rather than the absolute error is considered for convergence
+             * detection.
+             */
+            GmmxxNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
+
+            virtual AbstractNondeterministicLinearEquationSolver<Type>* clone() const;
             
-            virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const override {
-                // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-                std::unique_ptr<gmm::csr_matrix<Type>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
-                
-                // Create vector for result of multiplication, which is reduced to the result vector after
-                // each multiplication.
-                std::vector<Type> multiplyResult(A.getRowCount());
-                
-                // Now perform matrix-vector multiplication as long as we meet the bound of the formula.
-                for (uint_fast64_t i = 0; i < n; ++i) {
-                    gmm::mult(*gmmxxMatrix, x, multiplyResult);
-                    
-                    if (b != nullptr) {
-                        gmm::add(*b, multiplyResult);
-                    }
-                    
-                    if (minimize) {
-                        storm::utility::vector::reduceVectorMin(multiplyResult, x, nondeterministicChoiceIndices);
-                    } else {
-                        storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices);
-                    }
-                }
-            }
+            virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* b = nullptr, uint_fast64_t n = 1, std::vector<Type>* multiplyResult = nullptr) const override;
             
-            virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* multiplyResult = nullptr, std::vector<Type>* newX = nullptr) const override {
-                // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-                std::unique_ptr<gmm::csr_matrix<Type>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
-                
-                // Set up the environment for the power method.
-                bool multiplyResultMemoryProvided = true;
-                if (multiplyResult == nullptr) {
-                    multiplyResult = new std::vector<Type>(A.getRowCount());
-                    multiplyResultMemoryProvided = false;
-                }
-        
-                std::vector<Type>* currentX = &x;
-                bool xMemoryProvided = true;
-                if (newX == nullptr) {
-                    newX = new std::vector<Type>(x.size());
-                    xMemoryProvided = false;
-                }
-                std::vector<Type>* swap = nullptr;
-                uint_fast64_t iterations = 0;
-                bool converged = false;
-                
-                // Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number
-                // of iterations.
-                while (!converged && iterations < this->maxIterations) {
-                    // Compute x' = A*x + b.
-                    gmm::mult(*gmmxxMatrix, *currentX, *multiplyResult);
-                    gmm::add(b, *multiplyResult);
-                    
-                    // Reduce the vector x' by applying min/max for all non-deterministic choices.
-                    if (minimize) {
-                        storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, nondeterministicChoiceIndices);
-                    } else {
-                        storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, nondeterministicChoiceIndices);
-                    }
-                    
-                    // Determine whether the method converged.
-                    converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->precision, this->relative);
-                    
-                    // Update environment variables.
-                    swap = currentX;
-                    currentX = newX;
-                    newX = swap;
-                    ++iterations;
-                }
-                
-                // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
-                // is currently stored in currentX, but x is the output vector.
-                if (iterations % 2 == 1) {
-                    std::swap(x, *currentX);
-                    if (!xMemoryProvided) {
-                        delete currentX;
-                    }
-                } else if (!xMemoryProvided) {
-                    delete newX;
-                }
-        
-                if (!multiplyResultMemoryProvided) {
-                    delete multiplyResult;
-                }
-        
-                // Check if the solver converged and issue a warning otherwise.
-                if (converged) {
-                    LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
-                } else {
-                    LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations.");
-                }
-            }
+            virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* multiplyResult = nullptr, std::vector<Type>* newX = nullptr) const override;
+            
+        private:
+            // The required precision for the iterative methods.
+            double precision;
             
+            // Sets whether the relative or absolute error is to be considered for convergence detection.
+            bool relative;
+            
+            // The maximal number of iterations to do before iteration is aborted.
+            uint_fast64_t maximalNumberOfIterations;
         };
-        
     } // namespace solver
 } // namespace storm
 

src/solver/NativeLinearEquationSolver.cpp

@@ -0,0 +1,153 @@
+#include "src/solver/NativeLinearEquationSolver.h"
+#include "src/settings/Settings.h"
+#include "src/utility/vector.h"
+#include "src/exceptions/InvalidStateException.h"
+
+bool NativeLinearEquationSolverOptionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* instance) -> bool {
+    // Offer all available methods as a command line option.
+	std::vector<std::string> methods;
+	methods.push_back("jacobi");
+	instance->addOption(storm::settings::OptionBuilder("NativeLinearEquationSolver", "nativelin", "", "The method to be used for solving linear equation systems with the native engine. Available are: jacobi.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("jacobi").build()).build());
+    
+	instance->addOption(storm::settings::OptionBuilder("NativeLinearEquationSolver", "maxiter", "i", "The maximal number of iterations to perform before iterative solving is aborted.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(10000).build()).build());
+    
+	instance->addOption(storm::settings::OptionBuilder("NativeLinearEquationSolver", "precision", "", "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+    
+	instance->addOption(storm::settings::OptionBuilder("NativeLinearEquationSolver", "absolute", "", "Whether the relative or the absolute error is considered for deciding convergence.").build());
+    
+	return true;
+});
+
+namespace storm {
+    namespace solver {
+        
+        template<typename ValueType>
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : method(method), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ValueType>
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver() {
+            // Get the settings object to customize linear solving.
+            storm::settings::Settings* settings = storm::settings::Settings::getInstance();
+            
+            // Get appropriate settings.
+            maximalNumberOfIterations = settings->getOptionByLongName("maxiter").getArgument(0).getValueAsUnsignedInteger();
+            precision = settings->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
+            relative = settings->isSet("absolute");
+            
+            // Determine the method to be used.
+            std::string const& methodAsString = settings->getOptionByLongName("nativelin").getArgument(0).getValueAsString();
+            if (methodAsString == "jacobi") {
+                method = JACOBI;
+            }
+        }
+        
+        template<typename ValueType>
+        AbstractLinearEquationSolver<ValueType>* NativeLinearEquationSolver<ValueType>::clone() const {
+            return new NativeLinearEquationSolver<ValueType>(*this);
+        }
+        
+        template<typename ValueType>
+        void NativeLinearEquationSolver<ValueType>::solveEquationSystem(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+            // Get a Jacobi decomposition of the matrix A.
+            std::pair<storm::storage::SparseMatrix<ValueType>, storm::storage::SparseMatrix<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+            
+            // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
+            bool multiplyResultProvided = true;
+            std::vector<ValueType>* nextX = multiplyResult;
+            if (nextX == nullptr) {
+                nextX = new std::vector<ValueType>(x.size());
+                multiplyResultProvided = false;
+            }
+            std::vector<ValueType> const* copyX = nextX;
+            std::vector<ValueType>* currentX = &x;
+            
+            // Target vector for precision calculation.
+            std::vector<ValueType> tmpX(x.size());
+            
+            // Set up additional environment variables.
+            uint_fast64_t iterationCount = 0;
+            bool converged = false;
+            
+            while (!converged && iterationCount < maximalNumberOfIterations) {
+                // Compute D^-1 * (b - LU * x) and store result in nextX.
+                jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
+                storm::utility::vector::scaleVectorInPlace(tmpX, -storm::utility::constantOne<ValueType>());
+                storm::utility::vector::addVectorsInPlace(tmpX, b);
+                jacobiDecomposition.second.multiplyWithVector(tmpX, *nextX);
+                
+                // Swap the two pointers as a preparation for the next iteration.
+                std::swap(nextX, currentX);
+                
+                // Now check if the process already converged within our precision.
+                converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, precision, relative);
+                
+                // Increase iteration count so we can abort if convergence is too slow.
+                ++iterationCount;
+            }
+            
+            // If the last iteration did not write to the original x we have to swap the contents, because the
+            // output has to be written to the input parameter x.
+            if (currentX == copyX) {
+                std::swap(x, *currentX);
+            }
+            
+            // If the vector for the temporary multiplication result was not provided, we need to delete it.
+            if (!multiplyResultProvided) {
+                delete copyX;
+            }
+        }
+        
+        template<typename ValueType>
+        void NativeLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+            // Set up some temporary variables so that we can just swap pointers instead of copying the result after
+            // each iteration.
+            std::vector<ValueType>* swap = nullptr;
+            std::vector<ValueType>* currentX = &x;
+            
+            bool multiplyResultProvided = true;
+            std::vector<ValueType>* nextX = multiplyResult;
+            if (nextX == nullptr) {
+                nextX = new std::vector<ValueType>(x.size());
+                multiplyResultProvided = false;
+            }
+            std::vector<ValueType> const* copyX = nextX;
+            
+            // Now perform matrix-vector multiplication as long as we meet the bound.
+            for (uint_fast64_t i = 0; i < n; ++i) {
+                A.multiplyWithVector(*currentX, *nextX);
+                std::swap(nextX, currentX);
+                
+                // If requested, add an offset to the current result vector.
+                if (b != nullptr) {
+                    storm::utility::vector::addVectorsInPlace(*currentX, *b);
+                }
+            }
+            
+            // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x,
+            // because the output is supposed to be stored in the input vector x.
+            if (currentX == copyX) {
+                std::swap(x, *currentX);
+            }
+            
+            // If the vector for the temporary multiplication result was not provided, we need to delete it.
+            if (!multiplyResultProvided) {
+                delete copyX;
+            }
+        }
+
+
+        template<typename ValueType>
+        std::string NativeLinearEquationSolver<ValueType>::methodToString() const {
+            if (method == JACOBI) {
+                return "jacobi";
+            } else {
+                throw storm::exceptions::InvalidStateException() << "Illegal method '" << method << "' set in NativeLinearEquationSolver.";
+            }
+        }
+        
+        // Explicitly instantiate the linear equation solver.
+        template class NativeLinearEquationSolver<double>;
+    }
+}

src/solver/NativeLinearEquationSolver.h

@@ -0,0 +1,65 @@
+#ifndef STORM_SOLVER_NATIVELINEAREQUATIONSOLVER_H_
+#define STORM_SOLVER_NATIVELINEAREQUATIONSOLVER_H_
+
+#include "AbstractLinearEquationSolver.h"
+
+namespace storm {
+    namespace solver {
+        
+        /*!
+         * A class that uses StoRM's native matrix operations to implement the AbstractLinearEquationSolver interface.
+         */
+        template<typename ValueType>
+        class NativeLinearEquationSolver : public AbstractLinearEquationSolver<ValueType> {
+        public:
+            // An enumeration specifying the available solution methods.
+            enum SolutionMethod {
+                JACOBI
+            };
+            
+            /*!
+             * Constructs a linear equation solver with parameters being set according to the settings object.
+             */
+            NativeLinearEquationSolver();
+
+            /*!
+             * Constructs a linear equation solver with the given parameters.
+             *
+             * @param method The method to use for linear equation solving.
+             * @param precision The precision to use for convergence detection.
+             * @param maximalNumberOfIterations The maximal number of iterations do perform before iteration is aborted.
+             * @param relative If set, the relative error rather than the absolute error is considered for convergence
+             * detection.
+             */
+            NativeLinearEquationSolver(SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
+            
+            virtual AbstractLinearEquationSolver<ValueType>* clone() const override;
+            
+            virtual void solveEquationSystem(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            
+            virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+
+        private:
+            /*!
+             * Retrieves the string representation of the solution method associated with this solver.
+             *
+             * @return The string representation of the solution method associated with this solver.
+             */
+            std::string methodToString() const;
+            
+            // The method to use for solving linear equation systems.
+            SolutionMethod method;
+            
+            // The required precision for the iterative methods.
+            double precision;
+            
+            // Sets whether the relative or absolute error is to be considered for convergence detection.
+            bool relative;
+            
+            // The maximal number of iterations to do before iteration is aborted.
+            uint_fast64_t maximalNumberOfIterations;
+        };
+    }
+}
+
+#endif /* STORM_SOLVER_NATIVELINEAREQUATIONSOLVER_H_ */

src/solver/NativeNondeterministicLinearEquationSolver.cpp

@@ -0,0 +1,142 @@
+#include "src/solver/NativeNondeterministicLinearEquationSolver.h"
+#include "src/settings/Settings.h"
+#include "src/utility/vector.h"
+
+bool NativeNondeterministicLinearEquationSolverOptionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* instance) -> bool {
+	instance->addOption(storm::settings::OptionBuilder("NativeNondeterminsticLinearEquationSolver", "maxiter", "i", "The maximal number of iterations to perform before iterative solving is aborted.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(10000).build()).build());
+    
+	instance->addOption(storm::settings::OptionBuilder("NativeNondeterminsticLinearEquationSolver", "precision", "", "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+    
+	instance->addOption(storm::settings::OptionBuilder("NativeNondeterminsticLinearEquationSolver", "absolute", "", "Whether the relative or the absolute error is considered for deciding convergence.").build());
+    
+	return true;
+});
+
+namespace storm {
+    namespace solver {
+        
+        template<typename ValueType>
+        NativeNondeterministicLinearEquationSolver<ValueType>::NativeNondeterministicLinearEquationSolver() {
+            // Get the settings object to customize solving.
+            storm::settings::Settings* settings = storm::settings::Settings::getInstance();
+            
+            // Get appropriate settings.
+            maximalNumberOfIterations = settings->getOptionByLongName("maxiter").getArgument(0).getValueAsUnsignedInteger();
+            precision = settings->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
+            relative = settings->isSet("absolute");
+        }
+        
+        template<typename ValueType>
+        NativeNondeterministicLinearEquationSolver<ValueType>::NativeNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ValueType>
+        AbstractNondeterministicLinearEquationSolver<ValueType>* NativeNondeterministicLinearEquationSolver<ValueType>::clone() const {
+            return new NativeNondeterministicLinearEquationSolver<ValueType>(*this);
+        }
+        
+        template<typename ValueType>
+        void NativeNondeterministicLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+            
+            // Set up the environment for the power method. If scratch memory was not provided, we need to create it.
+            bool multiplyResultMemoryProvided = true;
+            if (multiplyResult == nullptr) {
+                multiplyResult = new std::vector<ValueType>(A.getRowCount());
+                multiplyResultMemoryProvided = false;
+            }
+            std::vector<ValueType>* currentX = &x;
+            bool xMemoryProvided = true;
+            if (newX == nullptr) {
+                newX = new std::vector<ValueType>(x.size());
+                xMemoryProvided = false;
+            }
+            std::vector<ValueType>* swap = nullptr;
+            uint_fast64_t iterations = 0;
+            bool converged = false;
+
+            // Keep track of which of the vectors for x is the auxiliary copy.
+            std::vector<ValueType>* copyX = newX;
+
+            // Proceed with the iterations as long as the method did not converge or reach the
+            // user-specified maximum number of iterations.
+            while (!converged && iterations < maximalNumberOfIterations) {
+                // Compute x' = A*x + b.
+                A.multiplyWithVector(*currentX, *multiplyResult);
+                storm::utility::vector::addVectorsInPlace(*multiplyResult, b);
+                
+                // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
+                // element of the min/max operator stack.
+                if (minimize) {
+                    storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, nondeterministicChoiceIndices);
+                } else {
+                    storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, nondeterministicChoiceIndices);
+                }
+                
+                // Determine whether the method converged.
+                converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative);
+                
+                // Update environment variables.
+                std::swap(currentX, newX);
+                ++iterations;
+            }
+            
+            // Check if the solver converged and issue a warning otherwise.
+            if (converged) {
+                LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
+            } else {
+                LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations.");
+            }
+    
+            // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
+            // is currently stored in currentX, but x is the output vector.
+            if (currentX == copyX) {
+                std::swap(x, *currentX);
+            }
+            
+            if (!xMemoryProvided) {
+                delete copyX;
+            }
+            
+            if (!multiplyResultMemoryProvided) {
+                delete multiplyResult;
+            }
+        }
+        
+        template<typename ValueType>
+        void NativeNondeterministicLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+            
+            // If scratch memory was not provided, we need to create it.
+            bool multiplyResultMemoryProvided = true;
+            if (multiplyResult == nullptr) {
+                multiplyResult = new std::vector<ValueType>(A.getRowCount());
+                multiplyResultMemoryProvided = false;
+            }
+
+            // Now perform matrix-vector multiplication as long as we meet the bound of the formula.
+            for (uint_fast64_t i = 0; i < n; ++i) {
+                A.multiplyWithVector(x, *multiplyResult);
+                
+                // Add b if it is non-null.
+                if (b != nullptr) {
+                    storm::utility::vector::addVectorsInPlace(*multiplyResult, *b);
+                }
+                
+                // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
+                // element of the min/max operator stack.
+                if (minimize) {
+                    storm::utility::vector::reduceVectorMin(*multiplyResult, x, nondeterministicChoiceIndices);
+                } else {
+                    storm::utility::vector::reduceVectorMax(*multiplyResult, x, nondeterministicChoiceIndices);
+                }
+            }
+            
+            if (!multiplyResultMemoryProvided) {
+                delete multiplyResult;
+            }
+        }
+        
+        // Explicitly instantiate the solver.
+        template class NativeNondeterministicLinearEquationSolver<double>;
+    } // namespace solver
+} // namespace storm

src/solver/NativeNondeterministicLinearEquationSolver.h

@@ -0,0 +1,50 @@
+#ifndef STORM_SOLVER_NATIVENONDETERMINISTICLINEAREQUATIONSOLVER_H_
+#define STORM_SOLVER_NATIVENONDETERMINISTICLINEAREQUATIONSOLVER_H_
+
+#include "src/solver/AbstractNondeterministicLinearEquationSolver.h"
+
+namespace storm {
+    namespace solver {
+        
+        /*!
+         * A class that uses the gmm++ library to implement the AbstractNondeterminsticLinearEquationSolver interface.
+         */
+        template<class ValueType>
+        class NativeNondeterministicLinearEquationSolver : public AbstractNondeterministicLinearEquationSolver<ValueType> {
+        public:
+            /*!
+             * Constructs a nondeterministic linear equation solver with parameters being set according to the settings
+             * object.
+             */
+            NativeNondeterministicLinearEquationSolver();
+            
+            /*!
+             * Constructs a nondeterminstic linear equation solver with the given parameters.
+             *
+             * @param precision The precision to use for convergence detection.
+             * @param maximalNumberOfIterations The maximal number of iterations do perform before iteration is aborted.
+             * @param relative If set, the relative error rather than the absolute error is considered for convergence
+             * detection.
+             */
+            NativeNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
+            
+            virtual AbstractNondeterministicLinearEquationSolver<ValueType>* clone() const override;
+            
+            virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* newX = nullptr) const override;
+            
+            virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
+
+        private:
+            // The required precision for the iterative methods.
+            double precision;
+            
+            // Sets whether the relative or absolute error is to be considered for convergence detection.
+            bool relative;
+            
+            // The maximal number of iterations to do before iteration is aborted.
+            uint_fast64_t maximalNumberOfIterations;
+        };
+    } // namespace solver
+} // namespace storm
+
+#endif /* STORM_SOLVER_NATIVENONDETERMINISTICLINEAREQUATIONSOLVER_H_ */

src/storm.cpp

@@ -40,6 +40,7 @@
 #include "src/parser/PrctlParser.h"
 #include "src/utility/ErrorHandling.h"
 #include "src/formula/Prctl.h"
+#include "src/utility/vector.h"
 
 #include "src/settings/Settings.h"
 // Registers all standard options
@@ -462,7 +463,7 @@ int main(const int argc, const char* argv[]) {
                     LOG4CPLUS_INFO(logger, "Model is a Markov automaton.");
                     std::shared_ptr<storm::models::MarkovAutomaton<double>> markovAutomaton = parser.getModel<storm::models::MarkovAutomaton<double>>();
                     markovAutomaton->close();
-                    storm::modelchecker::csl::SparseMarkovAutomatonCslModelChecker<double> mc(*markovAutomaton, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::AbstractNondeterministicLinearEquationSolver<double>()));
+                    storm::modelchecker::csl::SparseMarkovAutomatonCslModelChecker<double> mc(*markovAutomaton);
 //                    std::cout << mc.checkExpectedTime(true, markovAutomaton->getLabeledStates("goal")) << std::endl;
 //                    std::cout << mc.checkExpectedTime(false, markovAutomaton->getLabeledStates("goal")) << std::endl;
 //                    std::cout << mc.checkLongRunAverage(true, markovAutomaton->getLabeledStates("goal")) << std::endl;

src/utility/StormOptions.cpp

@@ -34,15 +34,20 @@ bool storm::utility::StormOptions::optionsRegistered = storm::settings::Settings
     
 	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "fixDeadlocks", "", "If the model contains deadlock states, setting this option will insert self-loops for these states.").build());
     
-	std::vector<std::string> matrixLibraries;
-	matrixLibraries.push_back("gmm++");
-	matrixLibraries.push_back("native");
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "matrixLibrary", "m", "Sets which matrix library is preferred for numerical operations.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("matrixLibraryName", "The name of a matrix library. Valid values are gmm++ and native.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(matrixLibraries)).setDefaultValueString("gmm++").build()).build());
+	std::vector<std::string> linearEquationSolver;
+	linearEquationSolver.push_back("gmm++");
+	linearEquationSolver.push_back("native");
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "linsolver", "", "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++ and native.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(linearEquationSolver)).setDefaultValueString("gmm++").build()).build());
+
+    std::vector<std::string> nondeterministicLinearEquationSolver;
+	nondeterministicLinearEquationSolver.push_back("gmm++");
+	nondeterministicLinearEquationSolver.push_back("native");
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "linsolver", "", "Sets which solver is preferred for solving systems of linear equations arising from nondeterministic systems.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the solver to prefer. Available are: gmm++ and native.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(nondeterministicLinearEquationSolver)).setDefaultValueString("native").build()).build());
 
     std::vector<std::string> lpSolvers;
 	lpSolvers.push_back("gurobi");
 	lpSolvers.push_back("glpk");
-	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "lpSolver", "", "Sets which LP solver is preferred.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("LP solver name", "The name of an available LP solver. Valid values are gurobi and glpk.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(lpSolvers)).setDefaultValueString("glpk").build()).build());
+	settings->addOption(storm::settings::OptionBuilder("StoRM Main", "lpsolver", "", "Sets which LP solver is preferred.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("LP solver name", "The name of an available LP solver. Valid values are gurobi and glpk.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(lpSolvers)).setDefaultValueString("glpk").build()).build());
     
 	return true;
 });

src/utility/solver.cpp

@@ -1,7 +1,13 @@
 #include "src/utility/solver.h"
 
+#include "src/settings/Settings.h"
+
+#include "src/solver/NativeLinearEquationSolver.h"
 #include "src/solver/GmmxxLinearEquationSolver.h"
+
+#include "src/solver/NativeNondeterministicLinearEquationSolver.h"
 #include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
+
 #include "src/solver/GurobiLpSolver.h"
 #include "src/solver/GlpkLpSolver.h"
 
@@ -9,21 +15,23 @@ namespace storm {
     namespace utility {
         namespace solver {
             std::shared_ptr<storm::solver::LpSolver> getLpSolver(std::string const& name) {
-                std::string const& lpSolver = storm::settings::Settings::getInstance()->getOptionByLongName("lpSolver").getArgument(0).getValueAsString();
+                std::string const& lpSolver = storm::settings::Settings::getInstance()->getOptionByLongName("lpsolver").getArgument(0).getValueAsString();
                 if (lpSolver == "gurobi") {
                     return std::shared_ptr<storm::solver::LpSolver>(new storm::solver::GurobiLpSolver(name));
                 } else if (lpSolver == "glpk") {
                     return std::shared_ptr<storm::solver::LpSolver>(new storm::solver::GlpkLpSolver(name));
                 }
                 
-                throw storm::exceptions::InvalidSettingsException() << "No suitable LP-solver selected.";
+                throw storm::exceptions::InvalidSettingsException() << "No suitable LP solver selected.";
             }
             
             template<typename ValueType>
             std::shared_ptr<storm::solver::AbstractLinearEquationSolver<ValueType>> getLinearEquationSolver() {
-                std::string const& matrixLibrary = storm::settings::Settings::getInstance()->getOptionByLongName("matrixLibrary").getArgument(0).getValueAsString();
-                if (matrixLibrary == "gmm++") {
+                std::string const& linearEquationSolver = storm::settings::Settings::getInstance()->getOptionByLongName("linsolver").getArgument(0).getValueAsString();
+                if (linearEquationSolver == "gmm++") {
                     return std::shared_ptr<storm::solver::AbstractLinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>());
+                } else if (linearEquationSolver == "native") {
+                    return std::shared_ptr<storm::solver::AbstractLinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>());
                 }
                 
                 throw storm::exceptions::InvalidSettingsException() << "No suitable linear equation solver selected.";
@@ -31,11 +39,11 @@ namespace storm {
             
             template<typename ValueType>
             std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>> getNondeterministicLinearEquationSolver() {
-                std::string const& matrixLibrary = storm::settings::Settings::getInstance()->getOptionByLongName("matrixLibrary").getArgument(0).getValueAsString();
-                if (matrixLibrary == "gmm++") {
+                std::string const& nondeterministicLinearEquationSolver = storm::settings::Settings::getInstance()->getOptionByLongName("ndsolver").getArgument(0).getValueAsString();
+                if (nondeterministicLinearEquationSolver == "gmm++") {
                     return std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<ValueType>());
-                } else if (matrixLibrary == "native") {
-                    return std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>>(new storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>());
+                } else if (nondeterministicLinearEquationSolver == "native") {
+                    return std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>>(new storm::solver::NativeNondeterministicLinearEquationSolver<ValueType>());
                 }
                 
                 throw storm::exceptions::InvalidSettingsException() << "No suitable nondeterministic linear equation solver selected.";

src/utility/vector.cpp

@@ -0,0 +1,16 @@
+#include "src/utility/vector.h"
+
+template<typename ValueType>
+std::ostream& operator<<(std::ostream& out, std::vector<ValueType> const& vector) {
+    out << "vector (" << vector.size() << ") [ ";
+    for (uint_fast64_t i = 0; i < vector.size() - 1; ++i) {
+        out << vector[i] << ", ";
+    }
+    out << vector.back();
+    out << " ]";
+    return out;
+}
+
+// Explicitly instantiate functions.
+template std::ostream& operator<<(std::ostream& out, std::vector<double> const& vector);
+template std::ostream& operator<<(std::ostream& out, std::vector<uint_fast64_t> const& vector);

src/utility/vector.h

@@ -15,6 +15,9 @@
 #include "log4cplus/loggingmacros.h"
 extern log4cplus::Logger logger;
 
+template<typename ValueType>
+std::ostream& operator<<(std::ostream& out, std::vector<ValueType> const& vector);
+
 namespace storm {
     namespace utility {
         namespace vector {
@@ -129,34 +132,80 @@ namespace storm {
             }
             
             /*!
-             * Adds the two given vectors and writes the result into the first operand.
+             * Applies the given operation pointwise on the two given vectors and writes the result into the first
+             * vector.
              *
-             * @param target The first summand and target vector.
-             * @param summand The second summand.
+             * @param target The first operand and target vector.
+             * @param secondOperand The second operand.
              */
             template<class T>
-            void addVectorsInPlace(std::vector<T>& target, std::vector<T> const& summand) {
+            void applyPointwiseInPlace(std::vector<T>& target, std::vector<T> const& secondOperand, std::function<T (T const&, T const&)> function) {
 #ifdef DEBUG
                 if (target.size() != summand.size()) {
-                    throw storm::exceptions::InvalidArgumentException() << "Invalid call to storm::utility::vector::addVectorsInPlace: operand lengths mismatch.";
+                    throw storm::exceptions::InvalidArgumentException() << "Invalid call to storm::utility::vector::applyPointwiseInPlace: operand lengths mismatch.";
                 }
 #endif
 #ifdef STORM_HAVE_INTELTBB
                 tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
                                   [&](tbb::blocked_range<uint_fast64_t> const& range) {
-                                        uint_fast64_t firstRow = range.begin();
-                                        uint_fast64_t endRow = range.end();
-                                      
-                                        typename std::vector<T>::iterator targetIt = target.begin() + firstRow;
-                                        for (typename std::vector<T>::const_iterator summandIt = summand.begin() + firstRow, summandIte = summand.begin() + endRow; summandIt != summandIte; ++summandIt, ++targetIt) {
-                                                *targetIt += *summandIt;
-                                        }
+                                      std::transform(target.begin() + range.begin(), target.begin() + range.end(), secondOperand.begin() + range.begin(), target.begin() + range.begin(), function);
                                   });
 #else
-                std::transform(target.begin(), target.end(), summand.begin(), target.begin(), std::plus<T>());
+                std::transform(target.begin(), target.end(), secondOperand.begin(), secondOperand.begin(), function);
 #endif
             }
             
+            /*!
+             * Applies the given function pointwise on the given vector.
+             *
+             * @param target The vector to which to apply the function.
+             * @param function The function to apply.
+             */
+            template<class T>
+            void applyPointwiseInPlace(std::vector<T>& target, std::function<T (T const&)> function) {
+#ifdef STORM_HAVE_INTELTBB
+                tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
+                                  [&](tbb::blocked_range<uint_fast64_t> const& range) {
+                                      std::transform(target.begin() + range.begin(), target.begin() + range.end(), target.begin() + range.begin(), function);
+                                  });
+#else
+                std::transform(target.begin(), target.end(), target.begin(), function);
+#endif
+            }
+            
+            /*!
+             * Adds the two given vectors and writes the result into the first operand.
+             *
+             * @param target The first summand and target vector.
+             * @param summand The second summand.
+             */
+            template<class T>
+            void addVectorsInPlace(std::vector<T>& target, std::vector<T> const& summand) {
+                applyPointwiseInPlace<T>(target, summand, std::plus<T>());
+            }
+            
+            /*!
+             * Subtracts the two given vectors and writes the result into the first operand.
+             *
+             * @param target The first summand and target vector.
+             * @param summand The second summand.
+             */
+            template<class T>
+            void subtractVectorsInPlace(std::vector<T>& target, std::vector<T> const& summand) {
+                applyPointwiseInPlace<T>(target, summand, std::minus<T>());
+            }
+            
+            /*!
+             * Subtracts the two given vectors and writes the result into the first operand.
+             *
+             * @param target The first summand and target vector.
+             * @param summand The second summand.
+             */
+            template<class T>
+            void scaleVectorInPlace(std::vector<T>& target, T const& factor) {
+                applyPointwiseInPlace<T>(target, [&] (T const& argument) { return argument * factor; });
+            }
+            
             /*!
              * Reduces the given source vector by selecting an element according to the given filter out of each row group.
              *
@@ -371,7 +420,6 @@ namespace storm {
                 
                 return subVector;
             }
-            
         } // namespace vector
     } // namespace utility
 } // namespace storm

test/functional/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp

@@ -1,231 +0,0 @@
-#include "gtest/gtest.h"
-#include "storm-config.h"
-
-#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
-#include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
-#include "src/settings/Settings.h"
-#include "src/parser/AutoParser.h"
-
-TEST(GmmxxMdpPrctlModelCheckerTest, Dice) {
-	storm::settings::Settings* s = storm::settings::Settings::getInstance();
-	storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", "", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.trans.rew");
-    
-	ASSERT_EQ(parser.getType(), storm::models::MDP);
-
-	std::shared_ptr<storm::models::Mdp<double>> mdp = parser.getModel<storm::models::Mdp<double>>();
-
-	ASSERT_EQ(mdp->getNumberOfStates(), 169ull);
-	ASSERT_EQ(mdp->getNumberOfTransitions(), 436ull);
-
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
-
-	storm::property::prctl::Ap<double>* apFormula = new storm::property::prctl::Ap<double>("two");
-	storm::property::prctl::Eventually<double>* eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	storm::property::prctl::ProbabilisticNoBoundOperator<double>* probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
-    
-	std::vector<double> result = mc.checkNoBoundOperator(*probFormula);
-    
-	ASSERT_LT(std::abs(result[0] - 0.0277777612209320068), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete probFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("two");
-	eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
-    
-	result = mc.checkNoBoundOperator(*probFormula);
-    
-	ASSERT_LT(std::abs(result[0] - 0.0277777612209320068), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete probFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("three");
-	eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
-    
-	result = mc.checkNoBoundOperator(*probFormula);
-    
-	ASSERT_LT(std::abs(result[0] - 0.0555555224418640136), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete probFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("three");
-	eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
-    
-	result = mc.checkNoBoundOperator(*probFormula);
-    
-	ASSERT_LT(std::abs(result[0] - 0.0555555224418640136), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete probFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("four");
-	eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
-    
-	result = mc.checkNoBoundOperator(*probFormula);
-    
-	ASSERT_LT(std::abs(result[0] - 0.083333283662796020508), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete probFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("four");
-	eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
-    
-	result = mc.checkNoBoundOperator(*probFormula);
-    
-	ASSERT_LT(std::abs(result[0] - 0.083333283662796020508), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete probFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("done");
-	storm::property::prctl::ReachabilityReward<double>* reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	storm::property::prctl::RewardNoBoundOperator<double>* rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, true);
-    
-	result = mc.checkNoBoundOperator(*rewardFormula);
-    
-	ASSERT_LT(std::abs(result[0] - 7.3333294987678527832), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete rewardFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("done");
-	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, false);
-    
-	result = mc.checkNoBoundOperator(*rewardFormula);;
-    
-	ASSERT_LT(std::abs(result[0] - 7.3333294987678527832), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete rewardFormula;
-    
-	storm::parser::AutoParser<double> stateRewardParser(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.state.rew", "");
-    
-	ASSERT_EQ(stateRewardParser.getType(), storm::models::MDP);
-    
-	std::shared_ptr<storm::models::Mdp<double>> stateRewardMdp = stateRewardParser.getModel<storm::models::Mdp<double>>();
-    
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
-    
-	apFormula = new storm::property::prctl::Ap<double>("done");
-	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, true);
-    
-	result = stateRewardModelChecker.checkNoBoundOperator(*rewardFormula);
-    
-	ASSERT_LT(std::abs(result[0] - 7.3333294987678527832), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete rewardFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("done");
-	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, false);
-    
-	result = stateRewardModelChecker.checkNoBoundOperator(*rewardFormula);
-    
-	ASSERT_LT(std::abs(result[0] - 7.3333294987678527832), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete rewardFormula;
-    
-	storm::parser::AutoParser<double> stateAndTransitionRewardParser(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.state.rew", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.trans.rew");
-    
-	ASSERT_EQ(stateAndTransitionRewardParser.getType(), storm::models::MDP);
-    
-	std::shared_ptr<storm::models::Mdp<double>> stateAndTransitionRewardMdp = stateAndTransitionRewardParser.getModel<storm::models::Mdp<double>>();
-    
-    storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
-
-	apFormula = new storm::property::prctl::Ap<double>("done");
-	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, true);
-    
-	result = stateAndTransitionRewardModelChecker.checkNoBoundOperator(*rewardFormula);
-    
-	ASSERT_LT(std::abs(result[0] - (2.0 * 7.3333294987678527832)), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete rewardFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("done");
-	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, false);
-    
-	result = stateAndTransitionRewardModelChecker.checkNoBoundOperator(*rewardFormula);
-    
-	ASSERT_LT(std::abs(result[0] - (2.0 * 7.3333294987678527832)), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete rewardFormula;
-}
-
-TEST(GmmxxMdpPrctlModelCheckerTest, AsynchronousLeader) {
-	storm::settings::Settings* s = storm::settings::Settings::getInstance();
-	storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/mdp/asynchronous_leader/leader4.tra", STORM_CPP_BASE_PATH "/examples/mdp/asynchronous_leader/leader4.lab", "", STORM_CPP_BASE_PATH "/examples/mdp/asynchronous_leader/leader4.trans.rew");
-
-	ASSERT_EQ(parser.getType(), storm::models::MDP);
-
-	std::shared_ptr<storm::models::Mdp<double>> mdp = parser.getModel<storm::models::Mdp<double>>();
-
-	ASSERT_EQ(mdp->getNumberOfStates(), 3172ull);
-	ASSERT_EQ(mdp->getNumberOfTransitions(), 7144ull);
-
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
-
-	storm::property::prctl::Ap<double>* apFormula = new storm::property::prctl::Ap<double>("elected");
-	storm::property::prctl::Eventually<double>* eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	storm::property::prctl::ProbabilisticNoBoundOperator<double>* probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
-
-	std::vector<double> result = mc.checkNoBoundOperator(*probFormula);
-
-	ASSERT_LT(std::abs(result[0] - 1.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete probFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
-
-	result = mc.checkNoBoundOperator(*probFormula);
-
-	ASSERT_LT(std::abs(result[0] - 1.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete probFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	storm::property::prctl::BoundedEventually<double>* boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 25);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, true);
-
-	result = mc.checkNoBoundOperator(*probFormula);
-
-	ASSERT_LT(std::abs(result[0] - 0.0625), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete probFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 25);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, false);
-
-	result = mc.checkNoBoundOperator(*probFormula);
-
-	ASSERT_LT(std::abs(result[0] - 0.0625), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete probFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	storm::property::prctl::ReachabilityReward<double>* reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	storm::property::prctl::RewardNoBoundOperator<double>* rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, true);
-
-	result = mc.checkNoBoundOperator(*rewardFormula);;
-
-	ASSERT_LT(std::abs(result[0] - 4.28568908480604982), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete rewardFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, false);
-
-	result = mc.checkNoBoundOperator(*rewardFormula);;
-
-	ASSERT_LT(std::abs(result[0] - 4.2856904354441400784), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete rewardFormula;
-}

test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp

@@ -1,6 +1,7 @@
 #include "gtest/gtest.h"
 #include "storm-config.h"
 
+#include "src/solver/NativeNondeterministicLinearEquationSolver.h"
 #include "src/settings/Settings.h"
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/parser/AutoParser.h"
@@ -16,7 +17,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
 	ASSERT_EQ(mdp->getNumberOfStates(), 169ull);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), 436ull);
     
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::AbstractNondeterministicLinearEquationSolver<double>()));
+	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
     
 	storm::property::prctl::Ap<double>* apFormula = new storm::property::prctl::Ap<double>("two");
 	storm::property::prctl::Eventually<double>* eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
@@ -84,8 +85,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	result = mc.checkNoBoundOperator(*rewardFormula);
     
-	ASSERT_LT(std::abs(result[0] - 7.3333294987678527832), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
+	ASSERT_LT(std::abs(result[0] - 7.3333323), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
     
 	apFormula = new storm::property::prctl::Ap<double>("done");
@@ -94,8 +94,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	result = mc.checkNoBoundOperator(*rewardFormula);;
     
-	ASSERT_LT(std::abs(result[0] - 7.3333294987678527832), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
+	ASSERT_LT(std::abs(result[0] - 7.3333323), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
     
 	storm::parser::AutoParser<double> stateRewardParser(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.state.rew", "");
@@ -104,7 +103,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	std::shared_ptr<storm::models::Mdp<double>> stateRewardMdp = stateRewardParser.getModel<storm::models::Mdp<double>>();
     
-    storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::AbstractNondeterministicLinearEquationSolver<double>()));
+    storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
 
 	apFormula = new storm::property::prctl::Ap<double>("done");
 	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
@@ -112,8 +111,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	result = stateRewardModelChecker.checkNoBoundOperator(*rewardFormula);
     
-	ASSERT_LT(std::abs(result[0] - 7.3333294987678527832), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
+	ASSERT_LT(std::abs(result[0] - 7.3333323), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
     
 	apFormula = new storm::property::prctl::Ap<double>("done");
@@ -122,8 +120,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	result = stateRewardModelChecker.checkNoBoundOperator(*rewardFormula);
     
-	ASSERT_LT(std::abs(result[0] - 7.3333294987678527832), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
+	ASSERT_LT(std::abs(result[0] - 7.3333323), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
     
 	storm::parser::AutoParser<double> stateAndTransitionRewardParser(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.state.rew", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.trans.rew");
@@ -132,7 +129,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	std::shared_ptr<storm::models::Mdp<double>> stateAndTransitionRewardMdp = stateAndTransitionRewardParser.getModel<storm::models::Mdp<double>>();
     
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::AbstractNondeterministicLinearEquationSolver<double>()));
+	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
     
 	apFormula = new storm::property::prctl::Ap<double>("done");
 	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
@@ -140,8 +137,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	result = stateAndTransitionRewardModelChecker.checkNoBoundOperator(*rewardFormula);
     
-	ASSERT_LT(std::abs(result[0] - (2 * 7.3333294987678527832)), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
+	ASSERT_LT(std::abs(result[0] - 14.66666611), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
     
 	apFormula = new storm::property::prctl::Ap<double>("done");
@@ -150,8 +146,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	result = stateAndTransitionRewardModelChecker.checkNoBoundOperator(*rewardFormula);
     
-	ASSERT_LT(std::abs(result[0] - (2 * 7.3333294987678527832)), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
+	ASSERT_LT(std::abs(result[0] - 14.66666611), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
 }
 
@@ -166,7 +161,7 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
 	ASSERT_EQ(mdp->getNumberOfStates(), 3172ull);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), 7144ull);
 
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::AbstractNondeterministicLinearEquationSolver<double>()));
+	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
 
 	storm::property::prctl::Ap<double>* apFormula = new storm::property::prctl::Ap<double>("elected");
 	storm::property::prctl::Eventually<double>* eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
@@ -214,8 +209,7 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
 
 	result = mc.checkNoBoundOperator(*rewardFormula);;
 
-	ASSERT_LT(std::abs(result[0] - 4.28568908480604982), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
+	ASSERT_LT(std::abs(result[0] - 4.285709145), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
 
 	apFormula = new storm::property::prctl::Ap<double>("elected");
@@ -224,7 +218,6 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
 
 	result = mc.checkNoBoundOperator(*rewardFormula);;
 
-	ASSERT_LT(std::abs(result[0] - 4.2856904354441400784), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
+	ASSERT_LT(std::abs(result[0] - 4.285709145), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
 }

test/functional/solver/GmmxxLinearEquationSolverTest.cpp

@@ -0,0 +1,224 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/solver/GmmxxLinearEquationSolver.h"
+#include "src/settings/Settings.h"
+
+TEST(GmmxxLinearEquationSolver, SolveWithStandardOptions) {
+    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::GmmxxLinearEquationSolver<double> solver);
+    
+    storm::solver::GmmxxLinearEquationSolver<double> solver;
+    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}
+
+TEST(GmmxxLinearEquationSolver, gmres) {
+    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::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::GMRES, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE, true, 50));
+    
+    storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::GMRES, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE, true, 50);
+    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}
+
+TEST(GmmxxLinearEquationSolver, qmr) {
+    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::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::QMR, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE));
+    
+    storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::QMR, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE, 50);
+    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}
+
+TEST(GmmxxLinearEquationSolver, bicgstab) {
+    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::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::BICGSTAB, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE));
+    
+    storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::BICGSTAB, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE);
+    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}
+
+TEST(GmmxxLinearEquationSolver, jacobi) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, 4));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 2));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, -1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, 1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, -5));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 2));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, 2));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, 4));
+
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(3);
+    std::vector<double> b = {11, -16, 1};
+    
+    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::JACOBI, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE));
+    
+    storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::JACOBI, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE);
+    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}
+
+TEST(GmmxxLinearEquationSolver, gmresilu) {
+    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::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::GMRES, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::ILU, true, 50));
+    
+    storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::GMRES, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE, true, 50);
+    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}
+
+TEST(GmmxxLinearEquationSolver, gmresdiag) {
+    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::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::GMRES, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::DIAGONAL, true, 50));
+    
+    storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::GMRES, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::NONE, true, 50);
+    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}
+
+TEST(GmmxxLinearEquationSolver, 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::GmmxxLinearEquationSolver<double> solver;
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(A, x, nullptr, 4));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}

test/functional/solver/NativeLinearEquationSolverTest.cpp

@@ -0,0 +1,33 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/solver/NativeLinearEquationSolver.h"
+#include "src/settings/Settings.h"
+
+TEST(NativeLinearEquationSolver, SolveWithStandardOptions) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, 4));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 2));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, -1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, 1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, -5));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 2));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, 2));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, 4));
+    
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(3);
+    std::vector<double> b = {11, -16, 1};
+    
+    ASSERT_NO_THROW(storm::solver::NativeLinearEquationSolver<double> solver);
+    
+    storm::solver::NativeLinearEquationSolver<double> solver;
+    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
+}

test/performance/modelchecker/GmmxxMdpPrctModelCheckerTest.cpp

@@ -1,212 +0,0 @@
-#include "gtest/gtest.h"
-#include "storm-config.h"
-
-#include "src/settings/Settings.h"
-#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
-#include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
-#include "src/parser/AutoParser.h"
-
-TEST(GmmxxMdpPrctlModelCheckerTest, AsynchronousLeader) {
-	storm::settings::Settings* s = storm::settings::Settings::getInstance();
-	storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/mdp/asynchronous_leader/leader7.tra", STORM_CPP_BASE_PATH "/examples/mdp/asynchronous_leader/leader7.lab", "", STORM_CPP_BASE_PATH "/examples/mdp/asynchronous_leader/leader7.trans.rew");
-
-	ASSERT_EQ(parser.getType(), storm::models::MDP);
-
-	std::shared_ptr<storm::models::Mdp<double>> mdp = parser.getModel<storm::models::Mdp<double>>();
-
-	ASSERT_EQ(mdp->getNumberOfStates(), 2095783ull);
-	ASSERT_EQ(mdp->getNumberOfTransitions(), 7714385ull);
-
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
-
-	storm::property::prctl::Ap<double>* apFormula = new storm::property::prctl::Ap<double>("elected");
-	storm::property::prctl::Eventually<double>* eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	storm::property::prctl::ProbabilisticNoBoundOperator<double>* probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
-
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F elected] on asynchronous_leader/leader7...");
-	std::vector<double> result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-
-	ASSERT_LT(std::abs(result[0] - 1.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete probFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
-
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F elected] on asynchronous_leader/leader7...");
-	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-
-	ASSERT_LT(std::abs(result[0] - 1.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete probFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	storm::property::prctl::BoundedEventually<double>* boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 25ull);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, true);
-
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F<=25 elected] on asynchronous_leader/leader7...");
-	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-
-	ASSERT_LT(std::abs(result[0] - 0.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete probFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 25ull);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, false);
-
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F<=25 elected] on asynchronous_leader/leader7...");
-	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-
-	ASSERT_LT(std::abs(result[0] - 0.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete probFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	storm::property::prctl::ReachabilityReward<double>* reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	storm::property::prctl::RewardNoBoundOperator<double>* rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, true);
-
-    LOG4CPLUS_WARN(logger, "Model Checking Rmin=? [F elected] on asynchronous_leader/leader7...");
-	result = mc.checkNoBoundOperator(*rewardFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-
-	ASSERT_LT(std::abs(result[0] - 6.172433512), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete rewardFormula;
-
-	apFormula = new storm::property::prctl::Ap<double>("elected");
-	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, false);
-
-    LOG4CPLUS_WARN(logger, "Model Checking Rmax=? [F elected] on asynchronous_leader/leader7...");
-	result = mc.checkNoBoundOperator(*rewardFormula);
-    LOG4CPLUS_WARN(logger, "Done");
-
-	ASSERT_LT(std::abs(result[0] - 6.1724344), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
-	delete rewardFormula;
-}
-
-TEST(GmmxxMdpPrctlModelCheckerTest, Consensus) {
-	storm::settings::Settings* s = storm::settings::Settings::getInstance();
-
-	storm::parser::AutoParser<double> parser(STORM_CPP_BASE_PATH "/examples/mdp/consensus/coin4_6.tra", STORM_CPP_BASE_PATH "/examples/mdp/consensus/coin4_6.lab", STORM_CPP_BASE_PATH "/examples/mdp/consensus/coin4_6.steps.state.rew", "");
-    
-	ASSERT_EQ(parser.getType(), storm::models::MDP);
-    
-	std::shared_ptr<storm::models::Mdp<double>> mdp = parser.getModel<storm::models::Mdp<double>>();
-    
-	ASSERT_EQ(mdp->getNumberOfStates(), 63616ull);
-	ASSERT_EQ(mdp->getNumberOfTransitions(), 213472ull);
-    
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
-    
-    storm::property::prctl::Ap<double>* apFormula = new storm::property::prctl::Ap<double>("finished");
-	storm::property::prctl::Eventually<double>* eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
-	storm::property::prctl::ProbabilisticNoBoundOperator<double>* probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
-    
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F finished] on consensus/coin4_6...");
-	std::vector<double> result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-        
-	ASSERT_LT(std::abs(result[31168] - 1.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-    delete probFormula;
-    
-    apFormula = new storm::property::prctl::Ap<double>("finished");
-    storm::property::prctl::Ap<double>* apFormula2 = new storm::property::prctl::Ap<double>("all_coins_equal_0");
-    storm::property::prctl::And<double>* andFormula = new storm::property::prctl::And<double>(apFormula, apFormula2);
-	eventuallyFormula = new storm::property::prctl::Eventually<double>(andFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
-    
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F finished & all_coins_equal_0] on consensus/coin4_6...");
-	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 0.43742828319177884388579), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-    delete probFormula;
-    
-    apFormula = new storm::property::prctl::Ap<double>("finished");
-    apFormula2 = new storm::property::prctl::Ap<double>("all_coins_equal_1");
-    andFormula = new storm::property::prctl::And<double>(apFormula, apFormula2);
-    eventuallyFormula = new storm::property::prctl::Eventually<double>(andFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
-    
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F finished & all_coins_equal_1] on consensus/coin4_6...");
-	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 0.52932863686144482340267813924583606421947479248047), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-    delete probFormula;
-    
-    apFormula = new storm::property::prctl::Ap<double>("finished");
-    apFormula2 = new storm::property::prctl::Ap<double>("agree");
-    storm::property::prctl::Not<double>* notFormula = new storm::property::prctl::Not<double>(apFormula2);
-    andFormula = new storm::property::prctl::And<double>(apFormula, notFormula);
-    eventuallyFormula = new storm::property::prctl::Eventually<double>(andFormula);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
-    
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F finished & !agree] on consensus/coin4_6...");
-	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 0.1041409700076474653673841), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-    delete probFormula;
-    
-    apFormula = new storm::property::prctl::Ap<double>("finished");
-	storm::property::prctl::BoundedEventually<double>* boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 50ull);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, true);
-    
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F<=50 finished] on consensus/coin4_6...");
-	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 0.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-    delete probFormula;
-    
-    apFormula = new storm::property::prctl::Ap<double>("finished");
-	boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 50ull);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, false);
-    
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F<=50 finished] on consensus/coin4_6...");
-	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 0.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-    delete probFormula;
-    
-    apFormula = new storm::property::prctl::Ap<double>("finished");
-	storm::property::prctl::ReachabilityReward<double>* reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	storm::property::prctl::RewardNoBoundOperator<double>* rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, true);
-    
-    LOG4CPLUS_WARN(logger, "Model Checking Rmin=? [F finished] on consensus/coin4_6...");
-	result = mc.checkNoBoundOperator(*rewardFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 1725.5933133943854045), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete rewardFormula;
-    
-	apFormula = new storm::property::prctl::Ap<double>("finished");
-	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
-	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, false);
-    
-    LOG4CPLUS_WARN(logger, "Model Checking Rmax=? [F finished] on consensus/coin4_6...");
-	result = mc.checkNoBoundOperator(*rewardFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 2183.1424220082612919213715), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
-	delete rewardFormula;
-
-}

test/performance/modelchecker/SparseMdpPrctlModelCheckerTest.cpp

@@ -3,7 +3,7 @@
 
 #include "src/settings/Settings.h"
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
-#include "src/solver/AbstractNondeterministicLinearEquationSolver.h"
+#include "src/solver/NativeNondeterministicLinearEquationSolver.h"
 #include "src/parser/AutoParser.h"
 
 TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
@@ -17,78 +17,60 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
 	ASSERT_EQ(mdp->getNumberOfStates(), 2095783ull);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), 7714385ull);
 
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::AbstractNondeterministicLinearEquationSolver<double>()));
+	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
 
 	storm::property::prctl::Ap<double>* apFormula = new storm::property::prctl::Ap<double>("elected");
 	storm::property::prctl::Eventually<double>* eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
 	storm::property::prctl::ProbabilisticNoBoundOperator<double>* probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
 
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F elected] on asynchronous_leader/leader7...");
 	std::vector<double> result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
 
 	ASSERT_LT(std::abs(result[0] - 1.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
 	delete probFormula;
 
 	apFormula = new storm::property::prctl::Ap<double>("elected");
 	eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
 	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
 
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F elected] on asynchronous_leader/leader7...");
 	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
 
 	ASSERT_LT(std::abs(result[0] - 1.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
 	delete probFormula;
 
 	apFormula = new storm::property::prctl::Ap<double>("elected");
 	storm::property::prctl::BoundedEventually<double>* boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 25);
 	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, true);
 
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F<=25 elected] on asynchronous_leader/leader7...");
 	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
 
 	ASSERT_LT(std::abs(result[0] - 0.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
 	delete probFormula;
 
 	apFormula = new storm::property::prctl::Ap<double>("elected");
 	boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 25);
 	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, false);
 
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F<=25 elected] on asynchronous_leader/leader7...");
 	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
 
 	ASSERT_LT(std::abs(result[0] - 0.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-
 	delete probFormula;
 
 	apFormula = new storm::property::prctl::Ap<double>("elected");
 	storm::property::prctl::ReachabilityReward<double>* reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
 	storm::property::prctl::RewardNoBoundOperator<double>* rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, true);
 
-    LOG4CPLUS_WARN(logger, "Model Checking Rmin=? [F elected] on asynchronous_leader/leader7...");
 	result = mc.checkNoBoundOperator(*rewardFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-
-	ASSERT_LT(std::abs(result[0] - 6.172433512), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 
+	ASSERT_LT(std::abs(result[0] - 6.172489569), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
 
 	apFormula = new storm::property::prctl::Ap<double>("elected");
 	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
 	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, false);
 
-    LOG4CPLUS_WARN(logger, "Model Checking Rmax=? [F elected] on asynchronous_leader/leader7...");
 	result = mc.checkNoBoundOperator(*rewardFormula);
-    LOG4CPLUS_WARN(logger, "Done");
-
-	ASSERT_LT(std::abs(result[0] - 6.1724344), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 
+	ASSERT_LT(std::abs(result[0] - 6.17248915), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
 }
 
@@ -106,18 +88,15 @@ TEST(SparseMdpPrctlModelCheckerTest, Consensus) {
 	ASSERT_EQ(mdp->getNumberOfStates(), 63616ull);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), 213472ull);
     
-	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<double>>(new storm::solver::AbstractNondeterministicLinearEquationSolver<double>()));
+	storm::modelchecker::prctl::SparseMdpPrctlModelChecker<double> mc(*mdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
     
     storm::property::prctl::Ap<double>* apFormula = new storm::property::prctl::Ap<double>("finished");
 	storm::property::prctl::Eventually<double>* eventuallyFormula = new storm::property::prctl::Eventually<double>(apFormula);
 	storm::property::prctl::ProbabilisticNoBoundOperator<double>* probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
     
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F finished] on consensus/coin4_6...");
 	std::vector<double> result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 1), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
     
+	ASSERT_LT(std::abs(result[31168] - 1.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
     delete probFormula;
     
     apFormula = new storm::property::prctl::Ap<double>("finished");
@@ -126,12 +105,9 @@ TEST(SparseMdpPrctlModelCheckerTest, Consensus) {
 	eventuallyFormula = new storm::property::prctl::Eventually<double>(andFormula);
 	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
     
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F finished & all_coins_equal_0] on consensus/coin4_6...");
 	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 0.43742828319177884388579), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
     
+	ASSERT_LT(std::abs(result[31168] - 0.4372725189), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
     delete probFormula;
     
     apFormula = new storm::property::prctl::Ap<double>("finished");
@@ -140,12 +116,9 @@ TEST(SparseMdpPrctlModelCheckerTest, Consensus) {
     eventuallyFormula = new storm::property::prctl::Eventually<double>(andFormula);
 	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
     
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F finished & all_coins_equal_1] on consensus/coin4_6...");
 	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 0.52932863686144482340267813924583606421947479248047), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
     
+	ASSERT_LT(std::abs(result[31168] - 0.5291510935), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
     delete probFormula;
     
     apFormula = new storm::property::prctl::Ap<double>("finished");
@@ -155,60 +128,45 @@ TEST(SparseMdpPrctlModelCheckerTest, Consensus) {
     eventuallyFormula = new storm::property::prctl::Eventually<double>(andFormula);
 	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
     
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F finished & !agree] on consensus/coin4_6...");
 	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 0.1041409700076474653673841), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
     
+	ASSERT_LT(std::abs(result[31168] - 0.1039268793), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
     delete probFormula;
     
     apFormula = new storm::property::prctl::Ap<double>("finished");
 	storm::property::prctl::BoundedEventually<double>* boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 50ull);
 	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, true);
     
-    LOG4CPLUS_WARN(logger, "Model Checking Pmin=? [F<=50 finished] on consensus/coin4_6...");
 	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
     
 	ASSERT_LT(std::abs(result[31168] - 0.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
     delete probFormula;
     
     apFormula = new storm::property::prctl::Ap<double>("finished");
 	boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 50ull);
 	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, false);
     
-    LOG4CPLUS_WARN(logger, "Model Checking Pmax=? [F<=50 finished] on consensus/coin4_6...");
 	result = mc.checkNoBoundOperator(*probFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
     
 	ASSERT_LT(std::abs(result[31168] - 0.0), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
     delete probFormula;
 
     apFormula = new storm::property::prctl::Ap<double>("finished");
 	storm::property::prctl::ReachabilityReward<double>* reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
 	storm::property::prctl::RewardNoBoundOperator<double>* rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, true);
     
-    LOG4CPLUS_WARN(logger, "Model Checking Rmin=? [F finished] on consensus/coin4_6...");
 	result = mc.checkNoBoundOperator(*rewardFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
-    
-	ASSERT_LT(std::abs(result[31168] - 1725.5933133943854045), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
     
+	ASSERT_LT(std::abs(result[31168] - 1727.998607), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
 	delete rewardFormula;
     
 	apFormula = new storm::property::prctl::Ap<double>("finished");
 	reachabilityRewardFormula = new storm::property::prctl::ReachabilityReward<double>(apFormula);
 	rewardFormula = new storm::property::prctl::RewardNoBoundOperator<double>(reachabilityRewardFormula, false);
     
-    LOG4CPLUS_WARN(logger, "Model Checking Rmax=? [F finished] on consensus/coin4_6...");
 	result = mc.checkNoBoundOperator(*rewardFormula);
-    LOG4CPLUS_WARN(logger, "Done.");
     
 	ASSERT_LT(std::abs(result[31168] - 2183.1424220082612919213715), s->getOptionByLongName("precision").getArgument(0).getValueAsDouble());
-    
 	delete rewardFormula;
 
 }

test/performance/storm-performance-tests.cpp

@@ -34,7 +34,7 @@ void setUpLogging() {
  * Creates an empty settings object as the standard instance of the Settings class.
  */
 void createEmptyOptions() {
-    const char* newArgv[] = {"storm-performance-tests", "--maxIterations", "20000"};
+    const char* newArgv[] = {"storm-performance-tests", "--maxiter", "20000"};
 	storm::settings::Settings* s = storm::settings::Settings::getInstance();
 	try {
 		storm::settings::Settings::parse(3, newArgv);