#include "src/solver/SymbolicLinearEquationSolver.h" #include "src/storage/dd/DdManager.h" #include "src/storage/dd/Add.h" #include "src/settings/SettingsManager.h" #include "src/settings/modules/NativeEquationSolverSettings.h" namespace storm { namespace solver { template SymbolicLinearEquationSolver::SymbolicLinearEquationSolver(storm::dd::Add const& A, storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : A(A), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) { // Intentionally left empty. } template SymbolicLinearEquationSolver::SymbolicLinearEquationSolver(storm::dd::Add const& A, storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) : A(A), allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) { // Get the settings object to customize solving. storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings(); // Get appropriate settings. maximalNumberOfIterations = settings.getMaximalIterationCount(); precision = settings.getPrecision(); relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative; } template storm::dd::Add SymbolicLinearEquationSolver::solveEquationSystem(storm::dd::Add const& x, storm::dd::Add const& b) const { // Start by computing the Jacobi decomposition of the matrix A. storm::dd::Add diagonal = x.getDdManager().template getAddOne(); for (auto const& pair : rowColumnMetaVariablePairs) { diagonal *= x.getDdManager().template getIdentity(pair.first).equals(x.getDdManager().template getIdentity(pair.second)).template toAdd(); diagonal *= x.getDdManager().getRange(pair.first).template toAdd() * x.getDdManager().getRange(pair.second).template toAdd(); } diagonal *= allRows.template toAdd(); storm::dd::Add lu = diagonal.ite(this->A.getDdManager().template getAddZero(), this->A); storm::dd::Add dinv = diagonal / (diagonal * this->A); // Set up additional environment variables. storm::dd::Add xCopy = x; uint_fast64_t iterationCount = 0; bool converged = false; while (!converged && iterationCount < maximalNumberOfIterations) { storm::dd::Add xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs); storm::dd::Add tmp = lu.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables); tmp = b - tmp; tmp = tmp.swapVariables(this->rowColumnMetaVariablePairs); tmp = dinv.multiplyMatrix(tmp, this->columnMetaVariables); // Now check if the process already converged within our precision. converged = xCopy.equalModuloPrecision(tmp, precision, relative); // If the method did not converge yet, we prepare the x vector for the next iteration. if (!converged) { xCopy = tmp; } // Increase iteration count so we can abort if convergence is too slow. ++iterationCount; } return xCopy; } template storm::dd::Add SymbolicLinearEquationSolver::performMatrixVectorMultiplication(storm::dd::Add const& x, storm::dd::Add const* b, uint_fast64_t n) const { storm::dd::Add xCopy = x; // Perform matrix-vector multiplication while the bound is met. for (uint_fast64_t i = 0; i < n; ++i) { xCopy = xCopy.swapVariables(this->rowColumnMetaVariablePairs); xCopy = this->A.multiplyMatrix(xCopy, this->columnMetaVariables); if (b != nullptr) { xCopy += *b; } } return xCopy; } template class SymbolicLinearEquationSolver; } }