diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index 123e0c440..e297d6506 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -55,7 +55,7 @@ namespace storm { } template - void GmmxxLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + bool GmmxxLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner '" << preconditionerToString() << "' (max. " << maximalNumberOfIterations << " iterations)."); if (method == SolutionMethod::Jacobi && preconditioner != Preconditioner::None) { LOG4CPLUS_WARN(logger, "Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored."); @@ -94,8 +94,10 @@ namespace storm { // Check if the solver converged and issue a warning otherwise. if (iter.converged()) { LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); + return true; } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + LOG4CPLUS_INFO(logger, "Iterative solver did not converge."); + return false; } } else if (method == SolutionMethod::Jacobi) { uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(*originalA, x, b, multiplyResult); @@ -103,9 +105,14 @@ namespace storm { // Check if the solver converged and issue a warning otherwise. if (iterations < maximalNumberOfIterations) { LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); + return true; } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + LOG4CPLUS_INFO(logger, "Iterative solver did not converge."); + return false; } + } else { + STORM_LOG_ERROR("The selected method is not supported."); + return false; } } diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index ec7aceb6d..bb446239c 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -46,10 +46,22 @@ namespace storm { */ GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A); - virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; + virtual bool solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; + virtual void setPrecision(double precision) override{ + this->precision = precision; + } + + virtual void setIterations(uint_fast64_t maximalNumberOfIterations) override{ + this->maximalNumberOfIterations = maximalNumberOfIterations; + } + + virtual void setRelative(bool relative) override{ + this->relative = relative; + } + private: /*! * Solves the linear equation system A*x = b given by the parameters using the Jacobi method. diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h index 61f26294e..d9dfafbe6 100644 --- a/src/solver/LinearEquationSolver.h +++ b/src/solver/LinearEquationSolver.h @@ -29,8 +29,10 @@ namespace storm { * @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A. * @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 true iff solving was successful (e.g. iterative methods converged) */ - virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const = 0; + virtual bool solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const = 0; /*! * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After @@ -45,6 +47,12 @@ namespace storm { * vector must be equal to the number of rows of A. */ virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const = 0; + + virtual void setPrecision(double precision) = 0; + virtual void setIterations(uint_fast64_t maximalNumberOfIterations) = 0; + + virtual void setRelative(bool relative) = 0; + }; } // namespace solver diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp index 09a24972c..ffd475ca9 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -28,7 +28,7 @@ namespace storm { } template - void NativeLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + bool NativeLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { if (method == NativeLinearEquationSolverSolutionMethod::SOR || method == NativeLinearEquationSolverSolutionMethod::GaussSeidel) { // Define the omega used for SOR. ValueType omega = method == NativeLinearEquationSolverSolutionMethod::SOR ? storm::settings::nativeEquationSolverSettings().getOmega() : storm::utility::one(); @@ -66,6 +66,9 @@ namespace storm { if (!tmpXProvided) { delete tmpX; } + + return converged; + } else { // Get a Jacobi decomposition of the matrix A. std::pair, std::vector> jacobiDecomposition = A.getJacobiDecomposition(); @@ -113,6 +116,8 @@ namespace storm { if (!multiplyResultProvided) { delete copyX; } + + return converged; } } diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index 6fa972f05..b177f7b3a 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -38,9 +38,22 @@ namespace storm { */ NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true); - virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; + virtual bool solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; + + virtual void setPrecision(double precision) override{ + this->precision = precision; + } + + virtual void setIterations(uint_fast64_t maximalNumberOfIterations) override{ + this->maximalNumberOfIterations = maximalNumberOfIterations; + } + + virtual void setRelative(bool relative) override{ + this->relative = relative; + } + private: /*! diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index 332641ceb..4ab6797b5 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -268,7 +268,6 @@ namespace storm { template SparseMatrix::SparseMatrix(SparseMatrix const& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(other.columnsAndValues), rowIndications(other.rowIndications), nontrivialRowGrouping(other.nontrivialRowGrouping), rowGroupIndices(other.rowGroupIndices) { // Intentionally left empty. - std::cout << "!!!! copying matrix (constructor)" << std::endl; } template @@ -276,7 +275,6 @@ namespace storm { storm::storage::BitVector rowConstraint(other.getRowCount(), true); storm::storage::BitVector columnConstraint(other.getColumnCount(), true); *this = other.getSubmatrix(false, rowConstraint, columnConstraint, insertDiagonalElements); - std::cout << "!!!! copying matrix (insertingDiagElements)" << std::endl; } template @@ -290,7 +288,6 @@ namespace storm { template SparseMatrix::SparseMatrix(index_type columnCount, std::vector const& rowIndications, std::vector> const& columnsAndValues, std::vector const& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(rowGroupIndices) { this->updateNonzeroEntryCount(); - std::cout << "!!!! copying matrix (copying Ingredients)" << std::endl; } template @@ -300,7 +297,6 @@ namespace storm { template SparseMatrix& SparseMatrix::operator=(SparseMatrix const& other) { - std::cout << "!!!! copying matrix (operator=)" << std::endl; // Only perform assignment if source and target are not the same. if (this != &other) { rowCount = other.rowCount; diff --git a/src/utility/policyguessing.cpp b/src/utility/policyguessing.cpp index f21e4577b..9c5b0c459 100644 --- a/src/utility/policyguessing.cpp +++ b/src/utility/policyguessing.cpp @@ -201,25 +201,30 @@ namespace storm { ){ //Get the submatrix/subvector A,x, and b and invoke linear equation solver storm::storage::SparseMatrix subA = A.getSubmatrix(true, probGreater0States, probGreater0States, true); - subA.convertToEquationSystem(); + storm::storage::SparseMatrix eqSysA(subA); + eqSysA.convertToEquationSystem(); std::vector subX(probGreater0States.getNumberOfSetBits()); storm::utility::vector::selectVectorValues(subX, probGreater0States, x); std::vector subB(probGreater0States.getNumberOfSetBits()); storm::utility::vector::selectVectorValues(subB, probGreater0States, b); - std::unique_ptr> linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory().create(subA); - linEqSysSolver->solveEquationSystem(subX, subB); + std::unique_ptr> linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory().create(eqSysA); + linEqSysSolver->setRelative(relative); - //Performing a couple of iterations makes the result "stable" when doing value iteration, i.e., - //if the given equation system is induced by an optimal policy, value iteration will terminate after the first iteration. - subA.convertToEquationSystem(); - linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory().create(subA); + std::size_t iterations = 0; std::vector copyX(subX.size()); - std::size_t iterations = 100000; //don't run into an endless loop... + ValueType newPrecision = precision; do { - linEqSysSolver->performMatrixVectorMultiplication(subX, &subB, 10, ©X); - --iterations; - } while(!storm::utility::vector::equalModuloPrecision(subX, copyX, precision, relative) && iterations>0); - STORM_LOG_WARN_COND(iterations>0, "Iterations on result of linear equation solver did not converge."); + linEqSysSolver->setPrecision(newPrecision); + if(!linEqSysSolver->solveEquationSystem(subX, subB)){ + break; //Solver did not converge.. so we have to go on with the current solution. + } + subA.multiplyWithVector(subX,copyX); + storm::utility::vector::addVectors(copyX, subB, copyX); // = Ax + b + ++iterations; + newPrecision *= 0.1; + } while(!storm::utility::vector::equalModuloPrecision(subX, copyX, precision, relative) && iterations<30); + + STORM_LOG_DEBUG("Required to increase the precision " << iterations << " times in order to obtain a precise result"); //fill in the result storm::utility::vector::setVectorValues(x, probGreater0States, subX); storm::utility::vector::setVectorValues(x, (~probGreater0States), prob0Value); diff --git a/src/utility/policyguessing.h b/src/utility/policyguessing.h index 0dd92d6b5..c5ec739ec 100644 --- a/src/utility/policyguessing.h +++ b/src/utility/policyguessing.h @@ -149,12 +149,17 @@ namespace storm { * It is not assumed that qualitative properties of the Graph defined by A have been checked, yet. * However, actual target states are already filtered out. * To ensure a unique solution, we also need to filter out the "prob0"-states. + * + * If the result does not satisfy "Ax+b = x (modulo precision)", the solver is executed + * again with increased precision. * * @param A the matrix of the equation system * @param x The initial guess of the solution. * @param b The vector of the equation system * @param targetChoices marks the rows in the matrix that have a positive probability to lead to a target state * @param prob0Value the value that is assigned to the states that have probability zero to reach a target + * @param const& precision The precision to be used by the solver + * @param relative sets whether to consider relative errors * @return The solution vector in the form of the vector x. */ template