diff --git a/src/storm/solver/GameSolver.cpp b/src/storm/solver/GameSolver.cpp index 599c6721b..d985fe5da 100644 --- a/src/storm/solver/GameSolver.cpp +++ b/src/storm/solver/GameSolver.cpp @@ -6,6 +6,12 @@ #include "storm/utility/vector.h" #include "storm/utility/graph.h" +#include "storm/solver/GmmxxLinearEquationSolver.h" +#include "storm/solver/EigenLinearEquationSolver.h" +#include "storm/solver/NativeLinearEquationSolver.h" +#include "storm/solver/EliminationLinearEquationSolver.h" + +#include "storm/exceptions/NotImplementedException.h" namespace storm { namespace solver { template @@ -37,16 +43,62 @@ namespace storm { auto const& pl1Row = player1Matrix.getRow(pl1State, this->player1SchedulerHint->getChoice(pl1State)); STORM_LOG_ASSERT(pl1Row.getNumberOfEntries()==1, "We assume that rows of player one have one entry."); uint_fast64_t pl2State = pl1Row.begin()->getColumn(); - selectedRows[pl1State] = player2Matrix.getRowGroupIndices()[pl2State] + this->player2SchedulerHint->getChoice(pl1State); + selectedRows[pl1State] = player2Matrix.getRowGroupIndices()[pl2State] + this->player2SchedulerHint->getChoice(pl2State); } //Get the matrix and the vector induced by this selection auto inducedMatrix = player2Matrix.selectRowsFromRowIndexSequence(selectedRows, true); inducedMatrix.convertToEquationSystem(); storm::utility::vector::selectVectorValues(tmpResult, selectedRows, b); + + // Solve the resulting equation system. + // Note that the linEqSolver might consider a slightly different interpretation of "equalModuloPrecision". Hence, we iteratively increase its precision. auto submatrixSolver = storm::solver::GeneralLinearEquationSolverFactory().create(std::move(inducedMatrix)); + submatrixSolver->setCachingEnabled(true); if (this->lowerBound) { submatrixSolver->setLowerBound(this->lowerBound.get()); } if (this->upperBound) { submatrixSolver->setUpperBound(this->upperBound.get()); } - submatrixSolver->solveEquations(x, tmpResult); + ValueType linEqPrecision = this->getPrecision(); + uint_fast64_t hintIterations = 0; + std::vector otherTmpResult(tmpResult.size()); + bool vectorsNotEqual = false; + do { + vectorsNotEqual = false; + // Set the precision of the underlying solver + auto* gmmSolver = dynamic_cast*>(submatrixSolver.get()); + auto* nativeSolver = dynamic_cast*>(submatrixSolver.get()); + auto* eigenSolver = dynamic_cast*>(submatrixSolver.get()); + auto* eliminationSolver = dynamic_cast*>(submatrixSolver.get()); + if (gmmSolver) { + auto newSettings = gmmSolver->getSettings(); + newSettings.setPrecision(linEqPrecision); + gmmSolver->setSettings(newSettings); + } else if (nativeSolver) { + auto newSettings = nativeSolver->getSettings(); + newSettings.setPrecision(linEqPrecision); + nativeSolver->setSettings(newSettings); + } else if (eigenSolver) { + eigenSolver->getSettings().setPrecision(linEqPrecision); + } else if (eliminationSolver) { + // elimination solver does not consider a precision so nothing to do + } else { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Changing the precision for the given linear equation solver type is not implemented"); + } + + // invoke the solver + submatrixSolver->solveEquations(x, tmpResult); + submatrixSolver->multiply(x, nullptr, otherTmpResult); + linEqPrecision *= storm::utility::convertNumber(0.5); + ++hintIterations; + auto otherTmpResIt = otherTmpResult.begin(); + auto xIt = x.begin(); + for(auto tmpResIt = tmpResult.begin(); tmpResIt != tmpResult.end(); ++tmpResIt) { + if (!storm::utility::vector::equalModuloPrecision(*tmpResIt + *xIt, *otherTmpResIt + *xIt, this->getPrecision(), this->getRelative())) { + vectorsNotEqual = true; + break; + } + ++otherTmpResIt; ++xIt; + } + } while (vectorsNotEqual); + STORM_LOG_WARN_COND(hintIterations == 1, "required " << hintIterations << " linear equation solver invokations to apply the scheduler hints in GameSolver"); } // Now perform the actual value iteration. diff --git a/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp b/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp index 0b665746a..e4ad81603 100644 --- a/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp @@ -12,6 +12,7 @@ #include "storm/utility/macros.h" #include "storm/exceptions/InvalidSettingsException.h" #include "storm/exceptions/InvalidStateException.h" +#include "storm/exceptions/NotImplementedException.h" namespace storm { namespace solver { @@ -112,6 +113,7 @@ namespace storm { // Create a solver that we will use throughout the procedure. We will modify the matrix in each iteration. auto solver = linearEquationSolverFactory->create(std::move(submatrix)); + setPrecisionRelativeOfLinearEquationSolver(solver, this->getPrecision(), this->getRelative()); if (this->lowerBound) { solver->setLowerBound(this->lowerBound.get()); } @@ -231,15 +233,42 @@ namespace storm { } if(this->schedulerHint) { - // Resolve the nondeterminism according to the scheduler hint and solve the resulting equation system. + // Resolve the nondeterminism according to the scheduler hint storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(this->schedulerHint->getChoices(), true); submatrix.convertToEquationSystem(); storm::utility::vector::selectVectorValues(*auxiliaryRowGroupVector, this->schedulerHint->getChoices(), this->A.getRowGroupIndices(), b); + // Solve the resulting equation system. + // Note that the linEqSolver might consider a slightly different interpretation of "equalModuloPrecision". Hence, we iteratively increase its precision. auto submatrixSolver = linearEquationSolverFactory->create(std::move(submatrix)); + submatrixSolver->setCachingEnabled(true); if (this->lowerBound) { submatrixSolver->setLowerBound(this->lowerBound.get()); } if (this->upperBound) { submatrixSolver->setUpperBound(this->upperBound.get()); } - submatrixSolver->solveEquations(x, *auxiliaryRowGroupVector); + ValueType linEqPrecision = this->getPrecision(); + // Temporarily shrink the auxiliary row vector + auxiliaryRowVector->resize(A.getRowGroupCount()); + auxiliaryRowVector->reserve(A.getRowCount()); + uint_fast64_t hintIterations = 0; + bool vectorsNotEqual = false; + do { + vectorsNotEqual = false; + setPrecisionRelativeOfLinearEquationSolver(submatrixSolver, linEqPrecision, this->getRelative()); + submatrixSolver->solveEquations(x, *auxiliaryRowGroupVector); + submatrixSolver->multiply(x, nullptr, *auxiliaryRowVector); + linEqPrecision *= storm::utility::convertNumber(0.5); + ++hintIterations; + auto otherTmpResIt = auxiliaryRowVector->begin(); + auto xIt = x.begin(); + for(auto tmpResIt = auxiliaryRowGroupVector->begin(); tmpResIt != auxiliaryRowGroupVector->end(); ++tmpResIt) { + if (!storm::utility::vector::equalModuloPrecision(*tmpResIt + *xIt, *otherTmpResIt + *xIt, this->getPrecision(), this->getRelative())) { + vectorsNotEqual = true; + break; + } + ++otherTmpResIt; ++xIt; + } + } while (vectorsNotEqual); + STORM_LOG_WARN_COND(hintIterations == 1, "required " << hintIterations << " linear equation solver invokations to apply the scheduler hint"); + auxiliaryRowVector->resize(A.getRowCount()); } std::vector* newX = auxiliaryRowGroupVector.get(); @@ -348,6 +377,38 @@ namespace storm { } } + + template + void StandardMinMaxLinearEquationSolver::setPrecisionRelativeOfLinearEquationSolver(std::unique_ptr>& linEqSolver, ValueType precision, bool relative) const { + auto* gmmSolver = dynamic_cast*>(linEqSolver.get()); + auto* nativeSolver = dynamic_cast*>(linEqSolver.get()); + auto* eigenSolver = dynamic_cast*>(linEqSolver.get()); + auto* eliminationSolver = dynamic_cast*>(linEqSolver.get()); + if (gmmSolver) { + auto newSettings = gmmSolver->getSettings(); + newSettings.setPrecision(precision); + newSettings.setRelativeTerminationCriterion(relative); + gmmSolver->setSettings(newSettings); + } else if (nativeSolver) { + auto newSettings = nativeSolver->getSettings(); + newSettings.setPrecision(precision); + newSettings.setRelativeTerminationCriterion(relative); + nativeSolver->setSettings(newSettings); + } else if (eigenSolver) { + eigenSolver->getSettings().setPrecision(precision); + // no relative flag for eigen solver + } else if (eliminationSolver) { + // elimination solver does not consider a precision so nothing to do + } else { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Changing the precision for the given linear equation solver type is not implemented"); + } + } + + template<> + void StandardMinMaxLinearEquationSolver::setPrecisionRelativeOfLinearEquationSolver(std::unique_ptr>& linEqSolver, storm::RationalNumber precision, bool relative) const { + // Intentionally left empty + } + template StandardMinMaxLinearEquationSolverSettings const& StandardMinMaxLinearEquationSolver::getSettings() const { return settings; diff --git a/src/storm/solver/StandardMinMaxLinearEquationSolver.h b/src/storm/solver/StandardMinMaxLinearEquationSolver.h index f7e1e40c6..0b33aef35 100644 --- a/src/storm/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/storm/solver/StandardMinMaxLinearEquationSolver.h @@ -66,6 +66,8 @@ namespace storm { Status updateStatusIfNotConverged(Status status, std::vector const& x, uint64_t iterations) const; void reportStatus(Status status, uint64_t iterations) const; + void setPrecisionRelativeOfLinearEquationSolver(std::unique_ptr>& linEqSolver, ValueType precision, bool relative) const; + /// The settings of this solver. StandardMinMaxLinearEquationSolverSettings settings; diff --git a/src/storm/storage/geometry/NativePolytope.cpp b/src/storm/storage/geometry/NativePolytope.cpp index 29f7b594f..6aceb06d5 100644 --- a/src/storm/storage/geometry/NativePolytope.cpp +++ b/src/storm/storage/geometry/NativePolytope.cpp @@ -3,7 +3,7 @@ #include "storm/utility/macros.h" #include "storm/utility/solver.h" -#include "storm/solver/Z3LPSolver.h" +#include "storm/solver/Z3LpSolver.h" #include "storm/solver/SmtSolver.h" #include "storm/storage/geometry/nativepolytopeconversion/QuickHull.h" #include "storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.h" @@ -42,7 +42,7 @@ namespace storm { } else { std::vector eigenPoints; eigenPoints.reserve(points.size()); - for(auto const& p : points){ + for (auto const& p : points){ eigenPoints.emplace_back(storm::adapters::EigenAdapter::toEigenVector(p)); }