Storm
A Modern Probabilistic Model Checker
Loading...
Searching...
No Matches
SymbolicNativeLinearEquationSolver.cpp
Go to the documentation of this file.
2
11#include "storm/utility/dd.h"
13
14namespace storm {
15namespace solver {
16
17template<storm::dd::DdType DdType, typename ValueType>
21
22template<storm::dd::DdType DdType, typename ValueType>
24 storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables,
25 std::set<storm::expressions::Variable> const& columnMetaVariables,
26 std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs)
27 : SymbolicNativeLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) {
28 this->setMatrix(A);
29}
30
31template<storm::dd::DdType DdType, typename ValueType>
33 storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables,
34 std::set<storm::expressions::Variable> const& columnMetaVariables,
35 std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs)
36 : SymbolicLinearEquationSolver<DdType, ValueType>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) {
37 // Intentionally left empty.
38}
39
40template<storm::dd::DdType DdType, typename ValueType>
41NativeLinearEquationSolverMethod SymbolicNativeLinearEquationSolver<DdType, ValueType>::getMethod(Environment const& env, bool isExactMode) const {
42 // Adjust the method if it is not supported by this solver.
43 // none was specified and we want exact computations
44 auto method = env.solver().native().getMethod();
45
46 if (isExactMode) {
47 if (method != NativeLinearEquationSolverMethod::RationalSearch) {
48 if (env.solver().native().isMethodSetFromDefault()) {
49 method = NativeLinearEquationSolverMethod::RationalSearch;
51 "Selecting '" + toString(method) +
52 "' as the solution technique to guarantee exact results. If you want to override this, please explicitly specify a different method.");
53 } else {
54 STORM_LOG_WARN("The selected solution method does not guarantee exact results.");
55 }
56 }
57 } else {
58 if (method != NativeLinearEquationSolverMethod::Power && method != NativeLinearEquationSolverMethod::RationalSearch &&
59 method != NativeLinearEquationSolverMethod::Jacobi) {
60 method = NativeLinearEquationSolverMethod::Jacobi;
61 STORM_LOG_INFO("The selected solution method is not supported in the dd engine. Falling back to '" + toString(method) + "'.");
62 }
63 STORM_LOG_WARN_COND_DEBUG(!env.solver().isForceSoundness(), "Sound computations are not supported in the Dd engine.");
64 }
65 return method;
66}
67
68template<storm::dd::DdType DdType, typename ValueType>
71 storm::dd::Add<DdType, ValueType> const& b) const {
72 switch (getMethod(env, std::is_same<ValueType, storm::RationalNumber>::value)) {
73 case NativeLinearEquationSolverMethod::Jacobi:
74 return solveEquationsJacobi(env, x, b);
75 case NativeLinearEquationSolverMethod::Power:
76 return solveEquationsPower(env, x, b);
77 case NativeLinearEquationSolverMethod::RationalSearch:
78 return solveEquationsRationalSearch(env, x, b);
79 default:
80 STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The selected solution technique is not supported.");
82 }
83}
84
85template<storm::dd::DdType DdType, typename ValueType>
88 storm::dd::DdManager<DdType>& manager = this->getDdManager();
89
90 ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision());
91 uint64_t maxIter = env.solver().native().getMaximalNumberOfIterations();
92 bool relative = env.solver().native().getRelativeTerminationCriterion();
93
94 STORM_LOG_INFO("Solving symbolic linear equation system with NativeLinearEquationSolver (jacobi)");
95
96 // Start by computing the Jacobi decomposition of the matrix A.
97 storm::dd::Bdd<DdType> diagonal = storm::utility::dd::getRowColumnDiagonal(x.getDdManager(), this->rowColumnMetaVariablePairs);
98 diagonal &= this->allRows;
99
100 storm::dd::Add<DdType, ValueType> lu = diagonal.ite(manager.template getAddZero<ValueType>(), this->A);
101 storm::dd::Add<DdType, ValueType> diagonalAdd = diagonal.template toAdd<ValueType>();
102 storm::dd::Add<DdType, ValueType> diag = diagonalAdd.multiplyMatrix(this->A, this->columnMetaVariables);
103
104 storm::dd::Add<DdType, ValueType> scaledLu = lu / diag;
105 storm::dd::Add<DdType, ValueType> scaledB = b / diag;
106
107 // Set up additional environment variables.
109 uint_fast64_t iterationCount = 0;
110 bool converged = false;
111
112 while (!converged && iterationCount < maxIter) {
113 storm::dd::Add<DdType, ValueType> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
114 storm::dd::Add<DdType, ValueType> tmp = scaledB - scaledLu.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables);
115
116 // Now check if the process already converged within our precision.
117 converged = tmp.equalModuloPrecision(xCopy, precision, relative);
118
119 xCopy = tmp;
120
121 // Increase iteration count so we can abort if convergence is too slow.
122 ++iterationCount;
123 }
124
125 if (converged) {
126 STORM_LOG_INFO("Iterative solver (jacobi) converged in " << iterationCount << " iterations.");
127 } else {
128 STORM_LOG_WARN("Iterative solver (jacobi) did not converge in " << iterationCount << " iterations.");
129 }
130
131 return xCopy;
132}
133
134template<storm::dd::DdType DdType, typename ValueType>
135typename SymbolicNativeLinearEquationSolver<DdType, ValueType>::PowerIterationResult
136SymbolicNativeLinearEquationSolver<DdType, ValueType>::performPowerIteration(storm::dd::Add<DdType, ValueType> const& x,
137 storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision,
138 bool relativeTerminationCriterion, uint64_t maximalIterations) const {
139 // Set up additional environment variables.
141 uint_fast64_t iterations = 0;
143
144 while (status == SolverStatus::InProgress && iterations < maximalIterations) {
145 storm::dd::Add<DdType, ValueType> currentXAsColumn = currentX.swapVariables(this->rowColumnMetaVariablePairs);
146 storm::dd::Add<DdType, ValueType> tmp = this->A.multiplyMatrix(currentXAsColumn, this->columnMetaVariables) + b;
147
148 // Now check if the process already converged within our precision.
149 if (tmp.equalModuloPrecision(currentX, precision, relativeTerminationCriterion)) {
151 }
152
153 // Set up next iteration.
154 ++iterations;
155 currentX = tmp;
156 }
157
158 return PowerIterationResult(status, iterations, currentX);
159}
160
161template<storm::dd::DdType DdType, typename ValueType>
162storm::dd::Add<DdType, ValueType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsPower(Environment const& env,
164 storm::dd::Add<DdType, ValueType> const& b) const {
165 STORM_LOG_INFO("Solving symbolic linear equation system with NativeLinearEquationSolver (power)");
166 ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision());
167 PowerIterationResult result =
168 performPowerIteration(x, b, precision, env.solver().native().getRelativeTerminationCriterion(), env.solver().native().getMaximalNumberOfIterations());
169
170 if (result.status == SolverStatus::Converged) {
171 STORM_LOG_INFO("Iterative solver (power iteration) converged in " << result.iterations << " iterations.");
172 } else {
173 STORM_LOG_WARN("Iterative solver (power iteration) did not converge in " << result.iterations << " iterations.");
174 }
175
176 return result.values;
177}
178
179template<storm::dd::DdType DdType, typename ValueType>
180bool SymbolicNativeLinearEquationSolver<DdType, ValueType>::isSolutionFixedPoint(storm::dd::Add<DdType, ValueType> const& x,
181 storm::dd::Add<DdType, ValueType> const& b) const {
182 storm::dd::Add<DdType, ValueType> xAsColumn = x.swapVariables(this->rowColumnMetaVariablePairs);
183 storm::dd::Add<DdType, ValueType> tmp = this->A.multiplyMatrix(xAsColumn, this->columnMetaVariables);
184 tmp += b;
185
186 return x == tmp;
187}
188
189template<storm::dd::DdType DdType, typename ValueType>
190template<typename RationalType, typename ImpreciseType>
191storm::dd::Add<DdType, RationalType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::sharpen(
192 uint64_t precision, SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, storm::dd::Add<DdType, ImpreciseType> const& x,
193 storm::dd::Add<DdType, RationalType> const& rationalB, bool& isSolution) {
195
196 for (uint64_t p = 1; p <= precision; ++p) {
197 sharpenedX = x.sharpenKwekMehlhorn(p);
198
199 isSolution = rationalSolver.isSolutionFixedPoint(sharpenedX, rationalB);
200 if (isSolution) {
201 break;
202 }
203 }
204
205 return sharpenedX;
206}
207
208template<storm::dd::DdType DdType, typename ValueType>
209template<typename RationalType, typename ImpreciseType>
210storm::dd::Add<DdType, RationalType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(
211 Environment const& env, SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver,
212 SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB,
214 // Storage for the rational sharpened vector and the power iteration intermediate vector.
217
218 // The actual rational search.
219 uint64_t overallIterations = 0;
220 uint64_t powerIterationInvocations = 0;
221 ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision());
222 uint64_t maxIter = env.solver().native().getMaximalNumberOfIterations();
223 bool relative = env.solver().native().getRelativeTerminationCriterion();
225 while (status == SolverStatus::InProgress && overallIterations < maxIter) {
226 typename SymbolicNativeLinearEquationSolver<DdType, ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(
227 currentX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter - overallIterations);
228
229 ++powerIterationInvocations;
230 STORM_LOG_TRACE("Completed " << powerIterationInvocations << " power iteration invocations, the last one with precision " << precision
231 << " completed in " << result.iterations << " iterations.");
232
233 // Count the iterations.
234 overallIterations += result.iterations;
235
236 // Compute maximal precision until which to sharpen.
237 uint64_t p =
238 storm::utility::convertNumber<uint64_t>(storm::utility::ceil(storm::utility::log10<ValueType>(storm::utility::one<ValueType>() / precision)));
239
240 bool isSolution = false;
241 sharpenedX = sharpen<RationalType, ImpreciseType>(p, rationalSolver, result.values, rationalB, isSolution);
242
243 if (isSolution) {
245 } else {
246 currentX = result.values;
247 precision /= storm::utility::convertNumber<ValueType, uint64_t>(10);
248 }
250 status = SolverStatus::Aborted;
251 }
252 }
253
254 if (status == SolverStatus::InProgress && overallIterations == maxIter) {
256 }
257
258 if (status == SolverStatus::Converged) {
259 STORM_LOG_INFO("Iterative solver (rational search) converged in " << overallIterations << " iterations.");
260 } else {
261 STORM_LOG_WARN("Iterative solver (rational search) did not converge in " << overallIterations << " iterations.");
262 }
263
264 return sharpenedX;
265}
266
267template<storm::dd::DdType DdType, typename ValueType>
268template<typename ImpreciseType>
269typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type
270SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(Environment const& env, storm::dd::Add<DdType, ValueType> const& x,
271 storm::dd::Add<DdType, ValueType> const& b) const {
272 return solveEquationsRationalSearchHelper<ValueType, ValueType>(env, *this, *this, b, this->getLowerBoundsVector(), b);
273}
274
275template<storm::dd::DdType DdType, typename ValueType>
276template<typename ImpreciseType>
277typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && !storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type
278SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(Environment const& env, storm::dd::Add<DdType, ValueType> const& x,
279 storm::dd::Add<DdType, ValueType> const& b) const {
280 storm::dd::Add<DdType, storm::RationalNumber> rationalB = b.template toValueType<storm::RationalNumber>();
281 SymbolicNativeLinearEquationSolver<DdType, storm::RationalNumber> rationalSolver(this->A.template toValueType<storm::RationalNumber>(), this->allRows,
282 this->rowMetaVariables, this->columnMetaVariables,
283 this->rowColumnMetaVariablePairs);
284
286 solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(env, rationalSolver, *this, rationalB, this->getLowerBoundsVector(), b);
287 return rationalResult.template toValueType<ValueType>();
288}
289
290template<storm::dd::DdType DdType, typename ValueType>
291template<typename ImpreciseType>
292typename std::enable_if<!std::is_same<ValueType, ImpreciseType>::value, storm::dd::Add<DdType, ValueType>>::type
293SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(Environment const& env, storm::dd::Add<DdType, ValueType> const& x,
294 storm::dd::Add<DdType, ValueType> const& b) const {
295 // First try to find a solution using the imprecise value type.
298 try {
299 impreciseX = this->getLowerBoundsVector().template toValueType<ImpreciseType>();
300 storm::dd::Add<DdType, ImpreciseType> impreciseB = b.template toValueType<ImpreciseType>();
301 SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> impreciseSolver(
302 this->A.template toValueType<ImpreciseType>(), this->allRows, this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs);
303
304 rationalResult = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(env, *this, impreciseSolver, b, impreciseX, impreciseB);
305 } catch (storm::exceptions::PrecisionExceededException const& e) {
306 STORM_LOG_WARN("Precision of value type was exceeded, trying to recover by switching to rational arithmetic.");
307
308 // Fall back to precise value type if the precision of the imprecise value type was exceeded.
309 rationalResult = solveEquationsRationalSearchHelper<ValueType, ValueType>(env, *this, *this, b, impreciseX.template toValueType<ValueType>(), b);
310 }
311 return rationalResult.template toValueType<ValueType>();
312}
313
314template<storm::dd::DdType DdType, typename ValueType>
315storm::dd::Add<DdType, ValueType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearch(
316 Environment const& env, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
317 STORM_LOG_INFO("Solving symbolic linear equation system with NativeLinearEquationSolver (rational search)");
318 return solveEquationsRationalSearchHelper<double>(env, x, b);
319}
320
321template<storm::dd::DdType DdType, typename ValueType>
323 auto method = getMethod(env, std::is_same<ValueType, storm::RationalNumber>::value);
324 if (method == NativeLinearEquationSolverMethod::Jacobi) {
326 }
328}
329
330template<storm::dd::DdType DdType, typename ValueType>
333 auto method = getMethod(env, std::is_same<ValueType, storm::RationalNumber>::value);
334 if (method == NativeLinearEquationSolverMethod::RationalSearch) {
335 requirements.requireLowerBounds();
336 }
337 return requirements;
338}
339
340template<storm::dd::DdType DdType, typename ValueType>
341std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicNativeLinearEquationSolverFactory<DdType, ValueType>::create(
342 Environment const& env) const {
343 return std::make_unique<SymbolicNativeLinearEquationSolver<DdType, ValueType>>();
344}
345
350
355} // namespace solver
356} // namespace storm
SolverEnvironment & solver()
storm::RationalNumber const & getPrecision() const
uint64_t const & getMaximalNumberOfIterations() const
bool const & getRelativeTerminationCriterion() const
storm::solver::NativeLinearEquationSolverMethod const & getMethod() const
NativeSolverEnvironment & native()
Add< LibraryType, ValueType > swapVariables(std::vector< std::pair< storm::expressions::Variable, storm::expressions::Variable > > const &metaVariablePairs) const
Swaps the given pairs of meta variables in the ADD.
Definition Add.cpp:292
Add< LibraryType, storm::RationalNumber > sharpenKwekMehlhorn(uint64_t precision) const
Retrieves the function that sharpens all values in the current ADD with the Kwek-Mehlhorn algorithm.
Definition Add.cpp:151
Add< LibraryType, ValueType > multiplyMatrix(Add< LibraryType, ValueType > const &otherMatrix, std::set< storm::expressions::Variable > const &summationMetaVariables) const
Multiplies the current ADD (representing a matrix) with the given matrix by summing over the given me...
Definition Add.cpp:372
bool equalModuloPrecision(Add< LibraryType, ValueType > const &other, ValueType const &precision, bool relative=true) const
Checks whether the current and the given ADD represent the same function modulo some given precision.
Definition Add.cpp:211
Bdd< LibraryType > ite(Bdd< LibraryType > const &thenBdd, Bdd< LibraryType > const &elseBdd) const
Performs an if-then-else with the given operands, i.e.
Definition Bdd.cpp:113
DdManager< LibraryType > & getDdManager() const
Retrieves the manager that is responsible for this DD.
Definition Dd.cpp:39
LinearEquationSolverRequirements & requireLowerBounds(bool critical=true)
An interface that represents an abstract symbolic linear equation solver.
void setMatrix(storm::dd::Add< DdType, ValueType > const &newA)
virtual std::unique_ptr< storm::solver::SymbolicLinearEquationSolver< DdType, ValueType > > create(Environment const &env) const override
An interface that represents an abstract symbolic linear equation solver.
virtual LinearEquationSolverRequirements getRequirements(Environment const &env) const override
Retrieves the requirements of the solver under the current settings.
virtual storm::dd::Add< DdType, ValueType > solveEquations(Environment const &env, storm::dd::Add< DdType, ValueType > const &x, storm::dd::Add< DdType, ValueType > const &b) const override
Solves the equation system A*x = b.
virtual LinearEquationSolverProblemFormat getEquationProblemFormat(Environment const &env) const override
Retrieves the format in which this solver expects to solve equations.
#define STORM_LOG_INFO(message)
Definition logging.h:29
#define STORM_LOG_WARN(message)
Definition logging.h:30
#define STORM_LOG_TRACE(message)
Definition logging.h:17
#define STORM_LOG_THROW(cond, exception, message)
Definition macros.h:30
#define STORM_LOG_WARN_COND_DEBUG(cond, message)
Definition macros.h:18
SFTBDDChecker::ValueType ValueType
std::string toString(GurobiSolverMethod const &method)
Yields a string representation of the GurobiSolverMethod.
storm::dd::Bdd< Type > getRowColumnDiagonal(storm::dd::DdManager< Type > const &ddManager, std::vector< std::pair< storm::expressions::Variable, storm::expressions::Variable > > const &rowColumnMetaVariablePairs)
Definition dd.cpp:91
bool isTerminate()
Check whether the program should terminate (due to some abort signal).
ValueType ceil(ValueType const &number)
LabParser.cpp.
Definition cli.cpp:18