diff --git a/src/storm/environment/solver/SolverEnvironment.cpp b/src/storm/environment/solver/SolverEnvironment.cpp index 3b53fced5..8371b8158 100644 --- a/src/storm/environment/solver/SolverEnvironment.cpp +++ b/src/storm/environment/solver/SolverEnvironment.cpp @@ -87,15 +87,30 @@ namespace storm { return linearEquationSolverTypeSetFromDefault; } + boost::optional SolverEnvironment::getPrecisionOfCurrentLinearEquationSolver() const { + switch (getLinearEquationSolverType()) { + case storm::solver::EquationSolverType::Gmmxx: + return gmmxx().getPrecision(); + case storm::solver::EquationSolverType::Eigen: + return eigen().getPrecision(); + case storm::solver::EquationSolverType::Native: + return native().getPrecision(); + case storm::solver::EquationSolverType::Elimination: + return boost::none; + default: + STORM_LOG_ASSERT(false, "The selected solver type is unknown."); + } + } + void SolverEnvironment::setLinearEquationSolverPrecision(storm::RationalNumber const& value) { STORM_LOG_ASSERT(getLinearEquationSolverType() == storm::solver::EquationSolverType::Native || getLinearEquationSolverType() == storm::solver::EquationSolverType::Gmmxx || getLinearEquationSolverType() == storm::solver::EquationSolverType::Eigen || getLinearEquationSolverType() == storm::solver::EquationSolverType::Elimination, "The current solver type is not respected in this method."); - native().setPrecision(value); - gmmxx().setPrecision(value); - eigen().setPrecision(value); + native().setPrecision(value); + gmmxx().setPrecision(value); + eigen().setPrecision(value); // Elimination solver does not have a precision } diff --git a/src/storm/environment/solver/SolverEnvironment.h b/src/storm/environment/solver/SolverEnvironment.h index beb94d428..ec591beb2 100644 --- a/src/storm/environment/solver/SolverEnvironment.h +++ b/src/storm/environment/solver/SolverEnvironment.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "storm/environment/Environment.h" #include "storm/environment/SubEnvironment.h" @@ -40,6 +41,7 @@ namespace storm { void setLinearEquationSolverType(storm::solver::EquationSolverType const& value); bool isLinearEquationSolverTypeSetFromDefaultValue() const; + boost::optional getPrecisionOfCurrentLinearEquationSolver() const; void setLinearEquationSolverPrecision(storm::RationalNumber const& value); void setLinearEquationSolverRelativeTerminationCriterion(bool value); diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index ddcfc64ca..ce7100571 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -5,6 +5,7 @@ #include "storm/environment/solver/MinMaxSolverEnvironment.h" #include "storm/utility/KwekMehlhorn.h" +#include "storm/utility/NumberTraits.h" #include "storm/utility/vector.h" #include "storm/utility/macros.h" @@ -52,7 +53,7 @@ namespace storm { template bool IterativeMinMaxLinearEquationSolver::internalSolveEquations(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { bool result = false; - switch (getMethod(env, std::is_same::value)) { + switch (getMethod(env, storm::NumberTraits::IsExact)) { case MinMaxMethod::ValueIteration: if (env.solver().isForceSoundness()) { result = solveEquationsSoundValueIteration(env, dir, x, b); @@ -112,13 +113,20 @@ namespace storm { // The solver that we will use throughout the procedure. std::unique_ptr> solver; + // The linear equation solver should be at least as precise as this solver + std::unique_ptr environmentOfSolver; + boost::optional precOfSolver = env.solver().getPrecisionOfCurrentLinearEquationSolver(); + if (!storm::NumberTraits::IsExact && precOfSolver && precOfSolver.get() > env.solver().minMax().getPrecision()) { + environmentOfSolver = std::make_unique(env); + environmentOfSolver->solver().setLinearEquationSolverPrecision(env.solver().minMax().getPrecision()); + } SolverStatus status = SolverStatus::InProgress; uint64_t iterations = 0; this->startMeasureProgress(); do { // Solve the equation system for the 'DTMC'. - solveInducedEquationSystem(env, solver, scheduler, x, subB, b); + solveInducedEquationSystem(environmentOfSolver ? *environmentOfSolver : env, solver, scheduler, x, subB, b); // Go through the multiplication result and see whether we can improve any of the choices. bool schedulerImproved = false; @@ -189,7 +197,7 @@ namespace storm { // Start by copying the requirements of the linear equation solver. MinMaxLinearEquationSolverRequirements requirements(this->linearEquationSolverFactory->getRequirements(env)); - auto method = getMethod(env, std::is_same::value); + auto method = getMethod(env, storm::NumberTraits::IsExact); if (method == MinMaxMethod::ValueIteration) { if (env.solver().isForceSoundness()) { // Interval iteration requires a unique solution and lower+upper bounds @@ -295,7 +303,15 @@ namespace storm { if (this->hasInitialScheduler()) { // Solve the equation system induced by the initial scheduler. std::unique_ptr> linEqSolver; - solveInducedEquationSystem(env, linEqSolver, this->getInitialScheduler(), x, *auxiliaryRowGroupVector, b); + // The linear equation solver should be at least as precise as this solver + std::unique_ptr environmentOfSolver; + boost::optional precOfSolver = env.solver().getPrecisionOfCurrentLinearEquationSolver(); + if (!storm::NumberTraits::IsExact && precOfSolver && precOfSolver.get() > env.solver().minMax().getPrecision()) { + environmentOfSolver = std::make_unique(env); + environmentOfSolver->solver().setLinearEquationSolverPrecision(env.solver().minMax().getPrecision()); + } + + solveInducedEquationSystem(environmentOfSolver ? *environmentOfSolver : env, linEqSolver, this->getInitialScheduler(), x, *auxiliaryRowGroupVector, b); // If we were given an initial scheduler and are maximizing (minimizing), our current solution becomes // always less-or-equal (greater-or-equal) than the actual solution. guarantee = maximize(dir) ? SolverGuarantee::LessOrEqual : SolverGuarantee::GreaterOrEqual; @@ -308,10 +324,10 @@ namespace storm { guarantee = SolverGuarantee::GreaterOrEqual; } } else if (this->hasCustomTerminationCondition()) { - if (this->getTerminationCondition().requiresGuarantee(SolverGuarantee::LessOrEqual)) { + if (this->getTerminationCondition().requiresGuarantee(SolverGuarantee::LessOrEqual) && this->hasLowerBound()) { this->createLowerBoundsVector(x); guarantee = SolverGuarantee::LessOrEqual; - } else if (this->getTerminationCondition().requiresGuarantee(SolverGuarantee::GreaterOrEqual)) { + } else if (this->getTerminationCondition().requiresGuarantee(SolverGuarantee::GreaterOrEqual) && this->hasUpperBound()) { this->createUpperBoundsVector(x); guarantee = SolverGuarantee::GreaterOrEqual; } diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 6f51c7d1c..6517daef0 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -4,6 +4,7 @@ #include "storm/utility/ConstantsComparator.h" #include "storm/utility/KwekMehlhorn.h" +#include "storm/utility/NumberTraits.h" #include "storm/utility/constants.h" #include "storm/utility/vector.h" #include "storm/exceptions/InvalidStateException.h" @@ -312,7 +313,7 @@ namespace storm { // Now check for termination. converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative); - terminate = this->terminateNow(*currentX, SolverGuarantee::LessOrEqual); + terminate = this->terminateNow(*currentX, guarantee); // Potentially show progress. this->showProgressIterative(iterations); @@ -332,7 +333,6 @@ namespace storm { template bool NativeLinearEquationSolver::solveEquationsPower(Environment const& env, std::vector& x, std::vector const& b) const { - STORM_LOG_THROW(this->hasLowerBound(), storm::exceptions::UnmetRequirementException, "Solver requires upper bound, but none was given."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (Power)"); // Prepare the solution vectors. @@ -340,12 +340,22 @@ namespace storm { this->cachedRowVector = std::make_unique>(getMatrixRowCount()); } std::vector* currentX = &x; + SolverGuarantee guarantee = SolverGuarantee::None; + if (this->hasCustomTerminationCondition()) { + if (this->getTerminationCondition().requiresGuarantee(SolverGuarantee::LessOrEqual) && this->hasLowerBound()) { + this->createLowerBoundsVector(*currentX); + guarantee = SolverGuarantee::LessOrEqual; + } else if (this->getTerminationCondition().requiresGuarantee(SolverGuarantee::GreaterOrEqual) && this->hasUpperBound()) { + this->createUpperBoundsVector(*currentX); + guarantee = SolverGuarantee::GreaterOrEqual; + } + } std::vector* newX = this->cachedRowVector.get(); // Forward call to power iteration implementation. this->startMeasureProgress(); ValueType precision = storm::utility::convertNumber(env.solver().native().getPrecision()); - PowerIterationResult result = this->performPowerIteration(currentX, newX, b, precision, env.solver().native().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, 0, env.solver().native().getMaximalNumberOfIterations(), env.solver().native().getPowerMethodMultiplicationStyle()); + PowerIterationResult result = this->performPowerIteration(currentX, newX, b, precision, env.solver().native().getRelativeTerminationCriterion(), guarantee, 0, env.solver().native().getMaximalNumberOfIterations(), env.solver().native().getPowerMethodMultiplicationStyle()); // Swap the result in place. if (currentX == this->cachedRowVector.get()) { @@ -387,10 +397,9 @@ namespace storm { template bool NativeLinearEquationSolver::solveEquationsSoundPower(Environment const& env, std::vector& x, std::vector const& b) const { - STORM_LOG_THROW(this->hasLowerBound(), storm::exceptions::UnmetRequirementException, "Solver requires upper bound, but none was given."); + STORM_LOG_THROW(this->hasLowerBound(), storm::exceptions::UnmetRequirementException, "Solver requires lower bound, but none was given."); STORM_LOG_THROW(this->hasUpperBound(), storm::exceptions::UnmetRequirementException, "Solver requires upper bound, but none was given."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (SoundPower)"); - std::vector* lowerX = &x; this->createLowerBoundsVector(*lowerX); @@ -833,7 +842,7 @@ namespace storm { template bool NativeLinearEquationSolver::internalSolveEquations(Environment const& env, std::vector& x, std::vector const& b) const { - switch(getMethod(env, std::is_same::value)) { + switch(getMethod(env, storm::NumberTraits::IsExact)) { case NativeLinearEquationSolverMethod::SOR: return this->solveEquationsSOR(env, x, b, storm::utility::convertNumber(env.solver().native().getSorOmega())); case NativeLinearEquationSolverMethod::GaussSeidel: @@ -911,7 +920,7 @@ namespace storm { template LinearEquationSolverProblemFormat NativeLinearEquationSolver::getEquationProblemFormat(Environment const& env) const { - auto method = getMethod(env, std::is_same::value); + auto method = getMethod(env, storm::NumberTraits::IsExact); if (method == NativeLinearEquationSolverMethod::Power || method == NativeLinearEquationSolverMethod::RationalSearch) { return LinearEquationSolverProblemFormat::FixedPointSystem; } else { @@ -922,7 +931,7 @@ namespace storm { template LinearEquationSolverRequirements NativeLinearEquationSolver::getRequirements(Environment const& env) const { LinearEquationSolverRequirements requirements; - auto method = getMethod(env, std::is_same::value); + auto method = getMethod(env, storm::NumberTraits::IsExact); if (method == NativeLinearEquationSolverMethod::Power && env.solver().isForceSoundness()) { requirements.requireBounds(); } else if (method == NativeLinearEquationSolverMethod::RationalSearch) { diff --git a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp index d9b7698de..037179e2c 100644 --- a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp @@ -241,7 +241,14 @@ namespace storm { } else { // If we were given an initial scheduler, we take its solution as the starting point. if (this->hasInitialScheduler()) { - localX = solveEquationsWithScheduler(env, this->getInitialScheduler(), x, b); + // The linear equation solver should be at least as precise as this solver + std::unique_ptr environmentOfSolver; + boost::optional precOfSolver = env.solver().getPrecisionOfCurrentLinearEquationSolver(); + if (!storm::NumberTraits::IsExact && precOfSolver && precOfSolver.get() > env.solver().minMax().getPrecision()) { + environmentOfSolver = std::make_unique(env); + environmentOfSolver->solver().setLinearEquationSolverPrecision(env.solver().minMax().getPrecision()); + } + localX = solveEquationsWithScheduler(environmentOfSolver ? *environmentOfSolver : env, this->getInitialScheduler(), x, b); } else { localX = this->getLowerBoundsVector(); } @@ -301,13 +308,21 @@ namespace storm { storm::dd::Bdd scheduler = this->hasInitialScheduler() ? this->getInitialScheduler() : this->A.sumAbstract(this->columnMetaVariables).maxAbstractRepresentative(this->choiceVariables); // Initialize linear equation solver. - std::unique_ptr> linearEquationSolver = linearEquationSolverFactory->create(env, this->allRows, this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs); + // It should be at least as precise as this solver. + std::unique_ptr environmentOfSolver; + boost::optional precOfSolver = env.solver().getPrecisionOfCurrentLinearEquationSolver(); + if (!storm::NumberTraits::IsExact && precOfSolver && precOfSolver.get() > env.solver().minMax().getPrecision()) { + environmentOfSolver = std::make_unique(env); + environmentOfSolver->solver().setLinearEquationSolverPrecision(env.solver().minMax().getPrecision()); + } + + std::unique_ptr> linearEquationSolver = linearEquationSolverFactory->create(environmentOfSolver ? *environmentOfSolver : env, this->allRows, this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs); this->forwardBounds(*linearEquationSolver); // Iteratively solve and improve the scheduler. uint64_t maxIter = env.solver().minMax().getMaximalNumberOfIterations(); while (!converged && iterations < maxIter) { - storm::dd::Add schedulerX = solveEquationsWithScheduler(env, *linearEquationSolver, scheduler, currentSolution, b, diagonal); + storm::dd::Add schedulerX = solveEquationsWithScheduler(environmentOfSolver ? *environmentOfSolver : env, *linearEquationSolver, scheduler, currentSolution, b, diagonal); // Policy improvement step. storm::dd::Add choiceValues = this->A.multiplyMatrix(schedulerX.swapVariables(this->rowColumnMetaVariablePairs), this->columnMetaVariables) + b;