diff --git a/resources/3rdparty/gmm-4.2/include/gmm/gmm_blas.h b/resources/3rdparty/gmm-4.2/include/gmm/gmm_blas.h index 1ec21bb10..032976c42 100644 --- a/resources/3rdparty/gmm-4.2/include/gmm/gmm_blas.h +++ b/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; } diff --git a/src/adapters/GmmxxAdapter.h b/src/adapters/GmmxxAdapter.h index 73897c639..624b43f80 100644 --- a/src/adapters/GmmxxAdapter.h +++ b/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); diff --git a/src/counterexamples/PathBasedSubsystemGenerator.h b/src/counterexamples/PathBasedSubsystemGenerator.h index c00fe0c4d..58678a390 100644 --- a/src/counterexamples/PathBasedSubsystemGenerator.h +++ b/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: diff --git a/src/solver/AbstractNondeterministicLinearEquationSolver.h b/src/solver/AbstractNondeterministicLinearEquationSolver.h index 1d2c93c10..08158742b 100644 --- a/src/solver/AbstractNondeterministicLinearEquationSolver.h +++ b/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 diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index dc8310f2b..402a33d59 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/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; } diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index 218fedf21..aaab4b31f 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/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. diff --git a/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp b/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp index 74d976171..f7c13526c 100644 --- a/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp +++ b/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; -}); \ No newline at end of file +}); + +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 diff --git a/src/solver/GmmxxNondeterministicLinearEquationSolver.h b/src/solver/GmmxxNondeterministicLinearEquationSolver.h index 33c5c97b6..54a48efb8 100644 --- a/src/solver/GmmxxNondeterministicLinearEquationSolver.h +++ b/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 diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp new file mode 100644 index 000000000..0176bf24f --- /dev/null +++ b/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>; + } +} \ No newline at end of file diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h new file mode 100644 index 000000000..ab3ce2d47 --- /dev/null +++ b/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_ */ diff --git a/src/solver/NativeNondeterministicLinearEquationSolver.cpp b/src/solver/NativeNondeterministicLinearEquationSolver.cpp new file mode 100644 index 000000000..df7edd137 --- /dev/null +++ b/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 diff --git a/src/solver/NativeNondeterministicLinearEquationSolver.h b/src/solver/NativeNondeterministicLinearEquationSolver.h new file mode 100644 index 000000000..4e5126254 --- /dev/null +++ b/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_ */ diff --git a/src/storm.cpp b/src/storm.cpp index c69418bc4..7600f3372 100644 --- a/src/storm.cpp +++ b/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; diff --git a/src/utility/StormOptions.cpp b/src/utility/StormOptions.cpp index 39e57a6b9..5d8eea24c 100644 --- a/src/utility/StormOptions.cpp +++ b/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; }); diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp index 762cc360b..93b811389 100644 --- a/src/utility/solver.cpp +++ b/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."; diff --git a/src/utility/vector.cpp b/src/utility/vector.cpp new file mode 100644 index 000000000..e4a1dd05c --- /dev/null +++ b/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); diff --git a/src/utility/vector.h b/src/utility/vector.h index a514efed2..f779f2f4a 100644 --- a/src/utility/vector.h +++ b/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 diff --git a/test/functional/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp deleted file mode 100644 index d7ac097b5..000000000 --- a/test/functional/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp +++ /dev/null @@ -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; -} diff --git a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp index 29b74594d..418044ac6 100644 --- a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp +++ b/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; } diff --git a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp new file mode 100644 index 000000000..077920f13 --- /dev/null +++ b/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()); +} \ No newline at end of file diff --git a/test/functional/solver/NativeLinearEquationSolverTest.cpp b/test/functional/solver/NativeLinearEquationSolverTest.cpp new file mode 100644 index 000000000..3ccfbbf08 --- /dev/null +++ b/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()); +} \ No newline at end of file diff --git a/test/performance/modelchecker/GmmxxMdpPrctModelCheckerTest.cpp b/test/performance/modelchecker/GmmxxMdpPrctModelCheckerTest.cpp deleted file mode 100644 index da34f7fec..000000000 --- a/test/performance/modelchecker/GmmxxMdpPrctModelCheckerTest.cpp +++ /dev/null @@ -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; - -} \ No newline at end of file diff --git a/test/performance/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/performance/modelchecker/SparseMdpPrctlModelCheckerTest.cpp index a9e489fd1..a95aba8d7 100644 --- a/test/performance/modelchecker/SparseMdpPrctlModelCheckerTest.cpp +++ b/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; } \ No newline at end of file diff --git a/test/performance/storm-performance-tests.cpp b/test/performance/storm-performance-tests.cpp index dc5152855..e06883d06 100644 --- a/test/performance/storm-performance-tests.cpp +++ b/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);