From 79730379e476b3f0cf2e7b9535d95fcd64ee11c7 Mon Sep 17 00:00:00 2001 From: dehnert Date: Sat, 14 Dec 2013 00:13:46 +0100 Subject: [PATCH] Started refactoring the linear equation system solvers. Former-commit-id: 72d647fd42c20fa5dc6975ff31d87eb45ce31936 --- src/adapters/GmmxxAdapter.h | 8 +- src/solver/AbstractLinearEquationSolver.h | 41 ++- ...ractNondeterministicLinearEquationSolver.h | 4 +- src/solver/GmmxxLinearEquationSolver.cpp | 268 ++++++++++++++- src/solver/GmmxxLinearEquationSolver.h | 322 +++++------------- ...xxNondeterministicLinearEquationSolver.cpp | 11 + ...mmxxNondeterministicLinearEquationSolver.h | 8 +- 7 files changed, 389 insertions(+), 273 deletions(-) create mode 100644 src/solver/GmmxxNondeterministicLinearEquationSolver.cpp diff --git a/src/adapters/GmmxxAdapter.h b/src/adapters/GmmxxAdapter.h index f0fb9cf9e..73897c639 100644 --- a/src/adapters/GmmxxAdapter.h +++ b/src/adapters/GmmxxAdapter.h @@ -31,12 +31,12 @@ public: * @return A pointer to a row-major sparse matrix in gmm++ format. */ template - static gmm::csr_matrix* toGmmxxSparseMatrix(storm::storage::SparseMatrix const& matrix) { + static std::unique_ptr> toGmmxxSparseMatrix(storm::storage::SparseMatrix const& matrix) { uint_fast64_t realNonZeros = matrix.getEntryCount(); LOG4CPLUS_DEBUG(logger, "Converting matrix with " << realNonZeros << " non-zeros to gmm++ format."); // Prepare the resulting matrix. - gmm::csr_matrix* result = new gmm::csr_matrix(matrix.rowCount, matrix.columnCount); + std::unique_ptr> result(new gmm::csr_matrix(matrix.getRowCount(), matrix.getColumnCount())); // Copy Row Indications std::copy(matrix.rowIndications.begin(), matrix.rowIndications.end(), result->jc.begin()); @@ -65,13 +65,13 @@ public: * @return A pointer to a row-major sparse matrix in gmm++ format. */ template - static gmm::csr_matrix* toGmmxxSparseMatrix(storm::storage::SparseMatrix&& matrix) { + static std::unique_ptr> toGmmxxSparseMatrix(storm::storage::SparseMatrix&& matrix) { uint_fast64_t realNonZeros = matrix.getEntryCount(); std::cout << "here?!" << std::endl; LOG4CPLUS_DEBUG(logger, "Converting matrix with " << realNonZeros << " non-zeros to gmm++ format."); // Prepare the resulting matrix. - gmm::csr_matrix* result = new gmm::csr_matrix(matrix.rowCount, matrix.columnCount); + std::unique_ptr> result(new gmm::csr_matrix(matrix.getRowCount(), matrix.getColumnCount())); typedef unsigned long long IND_TYPE; typedef std::vector vectorType_ull_t; diff --git a/src/solver/AbstractLinearEquationSolver.h b/src/solver/AbstractLinearEquationSolver.h index 25905e29e..39c5cab37 100644 --- a/src/solver/AbstractLinearEquationSolver.h +++ b/src/solver/AbstractLinearEquationSolver.h @@ -1,23 +1,52 @@ #ifndef STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_ #define STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_ -#include "src/storage/SparseMatrix.h" - #include +#include "src/storage/SparseMatrix.h" + namespace storm { namespace solver { + /*! + * A class that represents an abstract linear equation solver. In addition to solving a system of linear + * equations, the functionality to repeatedly multiply a matrix with a given vector is provided. + */ template class AbstractLinearEquationSolver { public: - + /*! + * Makes a copy of the linear equation solver. + * + * @return A pointer to a copy of the linear equation solver. + */ virtual AbstractLinearEquationSolver* clone() const = 0; - virtual void solveEquationSystem(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const = 0; - - virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1) const = 0; + /*! + * Solves the equation system A*x = b. The matrix A is required to be square and have a unique solution. + * The solution of the set of linear equations will be written to the vector x. + * + * @param A The coefficient matrix of the equation system. + * @param x The solution vector that has to be computed. Its length must be equal to the number of rows of A. + * @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. + */ + virtual void solveEquationSystem(storm::storage::SparseMatrix const& A, 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 + * performing the necessary multiplications, the result is written to the input vector x. + * + * @param A The matrix to use for matrix-vector multiplication. + * @param x The initial vector with which to perform matrix-vector multiplication. Its length must be equal + * to the number of rows of A. + * @param b If non-null, this vector is added after each multiplication. If given, 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. + */ + virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const = 0; }; } // namespace solver diff --git a/src/solver/AbstractNondeterministicLinearEquationSolver.h b/src/solver/AbstractNondeterministicLinearEquationSolver.h index 8d6b41acc..1d2c93c10 100644 --- a/src/solver/AbstractNondeterministicLinearEquationSolver.h +++ b/src/solver/AbstractNondeterministicLinearEquationSolver.h @@ -16,8 +16,8 @@ namespace storm { AbstractNondeterministicLinearEquationSolver() { storm::settings::Settings* s = storm::settings::Settings::getInstance(); precision = s->getOptionByLongName("precision").getArgument(0).getValueAsDouble(); - maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger(); - relative = s->getOptionByLongName("relative").getArgument(0).getValueAsBoolean(); + maxIterations = s->getOptionByLongName("maxiter").getArgument(0).getValueAsUnsignedInteger(); + relative = !s->isSet("absolute"); } AbstractNondeterministicLinearEquationSolver(double precision, uint_fast64_t maxIterations, bool relative) : precision(precision), maxIterations(maxIterations), relative(relative) { diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index f9173d735..dc8310f2b 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -2,27 +2,277 @@ #include #include +#include -bool GmmxxLinearEquationSolverOptionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* instance) -> bool { +#include "src/adapters/GmmxxAdapter.h" +#include "src/settings/Settings.h" +#include "src/utility/vector.h" +#include "src/utility/constants.h" +#include "src/exceptions/InvalidStateException.h" + +#include "gmm/gmm_matrix.h" +#include "gmm/gmm_iter_solvers.h" +bool GmmxxLinearEquationSolverOptionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* instance) -> bool { + // Offer all available methods as a command line option. std::vector methods; methods.push_back("bicgstab"); methods.push_back("qmr"); - methods.push_back("lscg"); methods.push_back("gmres"); methods.push_back("jacobi"); + instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "gmmlin", "", "The method to be used for solving linear equation systems with the gmm++ engine. Available are: bicgstab, qmr, gmres, jacobi.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("gmres").build()).build()); - std::vector preconditioner; + // Register available preconditioners. + std::vector preconditioner; preconditioner.push_back("ilu"); preconditioner.push_back("diagonal"); preconditioner.push_back("ildlt"); preconditioner.push_back("none"); + instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "gmmpre", "", "The preconditioning technique used for solving linear equation systems with the gmm++ engine. Available are: ilu, diagonal, none.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the preconditioning method.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(preconditioner)).setDefaultValueString("ilu").build()).build()); - instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "leMethod", "", "The Method used in Linear Equation Solving. Available are: bicgstab, qmr, lscg, gmres, jacobi").addArgument(storm::settings::ArgumentBuilder::createStringArgument("leMethodName", "The Name of the Method to use").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("gmres").build()).build()); - instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "preconditioner", "", "The Preconditioning Technique used in Linear Equation Solving. Available are: ilu, diagonal, ildlt, none").addArgument(storm::settings::ArgumentBuilder::createStringArgument("preconditionerName", "The Name of the Preconditioning Method").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(preconditioner)).setDefaultValueString("ilu").build()).build()); - instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "maxIterations", "", "Maximum number of Iterations to perform while solving a linear equation system").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("iterationCount", "Max. Iteration Count").setDefaultValueUnsignedInteger(10000).build()).build()); - instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "precision", "", "Precision used for iterative solving of linear equation systems").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("precisionValue", "Precision").setDefaultValueDouble(1e-6).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1000000.0)).build()).build()); - instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "relative", "", "Whether the relative or the absolute error is considered for deciding convergence via a given precision").addArgument(storm::settings::ArgumentBuilder::createBooleanArgument("useRelative", "relative or absolute comparison").setDefaultValueBoolean(true).build()).build()); + instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "gmmrestart", "", "The number of iteration until restarted methods are actually restarted.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The number of iterations.").setDefaultValueUnsignedInteger(50).build()).build()); + + instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "maxiter", "i", "The maximal number of iterations to perform before iterative solving is aborted.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(10000).build()).build()); + + instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "precision", "", "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build()); + + instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "absolute", "", "Whether the relative or the absolute error is considered for deciding convergence.").build()); return true; -}); \ No newline at end of file +}); + +namespace storm { + namespace solver { + + template + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(SolutionMethod method, double precision, bool relative, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, uint_fast64_t restart) : method(method), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner), restart(restart) { + // Intentionally left empty. + } + + + template + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver() { + // Get the settings object to customize linear solving. + storm::settings::Settings* settings = storm::settings::Settings::getInstance(); + + // Get appropriate settings. + maximalNumberOfIterations = settings->getOptionByLongName("maxiter").getArgument(0).getValueAsUnsignedInteger(); + precision = settings->getOptionByLongName("precision").getArgument(0).getValueAsDouble(); + relative = settings->isSet("absolute"); + restart = settings->getOptionByLongName("gmmrestart").getArgument(0).getValueAsUnsignedInteger(); + + // Determine the method to be used. + std::string const& methodAsString = settings->getOptionByLongName("gmmlin").getArgument(0).getValueAsString(); + if (methodAsString == "bicgstab") { + method = BICGSTAB; + } else if (methodAsString == "qmr") { + method = QMR; + } else if (methodAsString == "gmres") { + method = GMRES; + } else if (methodAsString == "jacobi") { + method = JACOBI; + } + + // Check which preconditioner to use. + std::string const& preconditionAsString = settings->getOptionByLongName("gmmpre").getArgument(0).getValueAsString(); + if (preconditionAsString == "ilu") { + preconditioner = ILU; + } else if (preconditionAsString == "diagonal") { + preconditioner = DIAGONAL; + } else if (preconditionAsString == "none") { + preconditioner = NONE; + } + } + + template + AbstractLinearEquationSolver* GmmxxLinearEquationSolver::clone() const { + return new GmmxxLinearEquationSolver(*this); + } + + template + void GmmxxLinearEquationSolver::solveEquationSystem(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner " << preconditionerToString() << "'."); + if (method == JACOBI && preconditioner != NONE) { + LOG4CPLUS_WARN(logger, "Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored."); + } + + if (method == BICGSTAB || method == QMR || method == GMRES) { + std::unique_ptr> gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + + // Prepare an iteration object that determines the accuracy and the maximum number of iterations. + gmm::iteration iter(precision, 0, maximalNumberOfIterations); + + if (method == BICGSTAB) { + if (preconditioner == ILU) { + gmm::bicgstab(*gmmA, x, b, gmm::ilu_precond>(*gmmA), iter); + } else if (preconditioner == DIAGONAL) { + gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), iter); + } else if (preconditioner == NONE) { + gmm::bicgstab(*gmmA, x, b, gmm::identity_matrix(), iter); + } + } else if (method == QMR) { + if (preconditioner == ILU) { + gmm::qmr(*gmmA, x, b, gmm::ilu_precond>(*gmmA), iter); + } else if (preconditioner == DIAGONAL) { + gmm::qmr(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), iter); + } else if (preconditioner == NONE) { + gmm::qmr(*gmmA, x, b, gmm::identity_matrix(), iter); + } + } else if (method == GMRES) { + if (preconditioner == ILU) { + gmm::gmres(*gmmA, x, b, gmm::ilu_precond>(*gmmA), restart, iter); + } else if (preconditioner == DIAGONAL) { + gmm::gmres(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), restart, iter); + } else if (preconditioner == NONE) { + gmm::gmres(*gmmA, x, b, gmm::identity_matrix(), restart, iter); + } + } + + // Check if the solver converged and issue a warning otherwise. + if (iter.converged()) { + LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + } + } else if (method == JACOBI) { + uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b, multiplyResult); + + // Check if the solver converged and issue a warning otherwise. + if (iterations < maximalNumberOfIterations) { + LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + } + } + } + + template + void GmmxxLinearEquationSolver::performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector* b, uint_fast64_t n, std::vector* multiplyResult) const { + // Transform the transition probability A to the gmm++ format to use its arithmetic. + std::unique_ptr> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + + // Set up some temporary variables so that we can just swap pointers instead of copying the result after + // each iteration. + std::vector* swap = nullptr; + std::vector* currentX = &x; + + bool multiplyResultProvided = true; + std::vector* nextX = multiplyResult; + if (nextX == nullptr) { + nextX = new std::vector(x.size()); + multiplyResultProvided = false; + } + std::vector const* copyX = nextX; + + // Now perform matrix-vector multiplication as long as we meet the bound. + for (uint_fast64_t i = 0; i < n; ++i) { + gmm::mult(*gmmxxMatrix, *currentX, *nextX); + std::swap(nextX, currentX); + + // If requested, add an offset to the current result vector. + if (b != nullptr) { + gmm::add(*b, *currentX); + } + } + + // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, + // because the output is supposed to be stored in the input vector x. + if (currentX == copyX) { + std::swap(x, *currentX); + } + + // If the vector for the temporary multiplication result was not provided, we need to delete it. + if (multiplyResultProvided) { + delete copyX; + } + } + + template + uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + // Get a Jacobi decomposition of the matrix A. + std::pair, storm::storage::SparseMatrix> jacobiDecomposition = A.getJacobiDecomposition(); + + // Convert the (inverted) diagonal matrix to gmm++'s format. + std::unique_ptr> gmmDinv = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.second)); + // Convert the LU matrix to gmm++'s format. + std::unique_ptr> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.first)); + + // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with. + bool multiplyResultProvided = true; + std::vector* nextX = multiplyResult; + if (nextX == nullptr) { + nextX = new std::vector(x.size()); + multiplyResultProvided = false; + } + std::vector const* copyX = nextX; + std::vector* currentX = &x; + + // Target vector for precision calculation. + std::vector tmpX(x.size()); + + // Set up additional environment variables. + uint_fast64_t iterationCount = 0; + bool converged = false; + + while (!converged && iterationCount < maximalNumberOfIterations) { + // Compute D^-1 * (b - LU * x) and store result in nextX. + gmm::mult(*gmmLU, *currentX, tmpX); + gmm::add(b, gmm::scaled(tmpX, -storm::utility::constantOne()), tmpX); + gmm::mult(*gmmDinv, tmpX, *nextX); + + // Swap the two pointers as a preparation for the next iteration. + std::swap(nextX, currentX); + + // Now check if the process already converged within our precision. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, precision, relative); + + // Increase iteration count so we can abort if convergence is too slow. + ++iterationCount; + } + + // If the last iteration did not write to the original x we have to swap the contents, because the + // output has to be written to the input parameter x. + if (currentX == copyX) { + std::swap(x, *currentX); + } + + // If the vector for the temporary multiplication result was not provided, we need to delete it. + if (!multiplyResultProvided) { + delete copyX; + } + + return iterationCount; + } + + template + std::string GmmxxLinearEquationSolver::methodToString() const { + if (method == BICGSTAB) { + return "bicgstab"; + } else if (method == QMR) { + return "qmr"; + } else if (method == GMRES) { + return "gmres"; + } else if (method == JACOBI) { + return "jacobi"; + } else { + throw storm::exceptions::InvalidStateException() << "Illegal method '" << method << "' set in GmmxxLinearEquationSolver."; + } + } + + template + std::string GmmxxLinearEquationSolver::preconditionerToString() const { + if (preconditioner == ILU) { + return "ilu"; + } else if (preconditioner == DIAGONAL) { + return "diagonal"; + } else if (preconditioner == NONE) { + return "none"; + } else { + throw storm::exceptions::InvalidStateException() << "Illegal preconditioner '" << preconditioner << "' set in GmmxxLinearEquationSolver."; + } + } + + // Explicitly instantiate the solver. + template class GmmxxLinearEquationSolver; + } +} \ No newline at end of file diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index 5f7e2d04e..218fedf21 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -2,270 +2,100 @@ #define STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_ #include "AbstractLinearEquationSolver.h" -#include "src/adapters/GmmxxAdapter.h" -#include "src/utility/constants.h" -#include "src/settings/Settings.h" -#include "src/utility/vector.h" - -#include "gmm/gmm_matrix.h" -#include "gmm/gmm_iter_solvers.h" - -#include namespace storm { namespace solver { - template - class GmmxxLinearEquationSolver : public AbstractLinearEquationSolver { + /*! + * A class that uses the gmm++ library to implement the AbstractLinearEquationSolver interface. + */ + template + class GmmxxLinearEquationSolver : public AbstractLinearEquationSolver { public: + // An enumeration specifying the available preconditioners. + enum Preconditioner { + ILU, DIAGONAL, NONE + }; - virtual AbstractLinearEquationSolver* clone() const { - return new GmmxxLinearEquationSolver(); - } + // An enumeration specifying the available solution methods. + enum SolutionMethod { + BICGSTAB, QMR, GMRES, JACOBI + }; - virtual void solveEquationSystem(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { - // Get the settings object to customize linear solving. - storm::settings::Settings* s = storm::settings::Settings::getInstance(); - - // Prepare an iteration object that determines the accuracy, maximum number of iterations - // and the like. - uint_fast64_t maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger(); - gmm::iteration iter(s->getOptionByLongName("precision").getArgument(0).getValueAsDouble(), 0, maxIterations); - - // Print some information about the used preconditioner. - std::string const precond = s->getOptionByLongName("preconditioner").getArgument(0).getValueAsString(); - LOG4CPLUS_INFO(logger, "Starting iterative solver."); - - // ALL available solvers must be declared in the cpp File, where the options are registered! - // Dito for the Preconditioners - - std::string const chosenLeMethod = s->getOptionByLongName("leMethod").getArgument(0).getValueAsString(); - if (chosenLeMethod == "jacobi") { - if (precond != "none") { - LOG4CPLUS_WARN(logger, "Requested preconditioner '" << precond << "', which is unavailable for the Jacobi method. Dropping preconditioner."); - } - } else { - if (precond == "ilu") { - LOG4CPLUS_INFO(logger, "Using ILU preconditioner."); - } else if (precond == "diagonal") { - LOG4CPLUS_INFO(logger, "Using diagonal preconditioner."); - } else if (precond == "ildlt") { - LOG4CPLUS_INFO(logger, "Using ILDLT preconditioner."); - } else if (precond == "none") { - LOG4CPLUS_INFO(logger, "Using no preconditioner."); - } - } - - // Now do the actual solving. - if (chosenLeMethod == "bicgstab") { - LOG4CPLUS_INFO(logger, "Using BiCGStab method."); - // Transform the transition probability matrix to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); - if (precond == "ilu") { - gmm::bicgstab(*gmmA, x, b, gmm::ilu_precond>(*gmmA), iter); - } else if (precond == "diagonal") { - gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), iter); - } else if (precond == "ildlt") { - gmm::bicgstab(*gmmA, x, b, gmm::ildlt_precond>(*gmmA), iter); - } else if (precond == "none") { - gmm::bicgstab(*gmmA, x, b, gmm::identity_matrix(), iter); - } - - // Check if the solver converged and issue a warning otherwise. - if (iter.converged()) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); - } - delete gmmA; - } else if (chosenLeMethod == "qmr") { - LOG4CPLUS_INFO(logger, "Using QMR method."); - // Transform the transition probability matrix to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); - if (precond == "ilu") { - gmm::qmr(*gmmA, x, b, gmm::ilu_precond>(*gmmA), iter); - } else if (precond == "diagonal") { - gmm::qmr(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), iter); - } else if (precond == "ildlt") { - gmm::qmr(*gmmA, x, b, gmm::ildlt_precond>(*gmmA), iter); - } else if (precond == "none") { - gmm::qmr(*gmmA, x, b, gmm::identity_matrix(), iter); - } - - // Check if the solver converged and issue a warning otherwise. - if (iter.converged()) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); - } - delete gmmA; - } else if (chosenLeMethod == "lscg") { - LOG4CPLUS_INFO(logger, "Using LSCG method."); - // Transform the transition probability matrix to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); - - if (precond != "none") { - LOG4CPLUS_WARN(logger, "Requested preconditioner '" << precond << "', which is unavailable for the LSCG method. Dropping preconditioner."); - } - gmm::least_squares_cg(*gmmA, x, b, iter); - - // Check if the solver converged and issue a warning otherwise. - if (iter.converged()) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); - } - delete gmmA; - } else if (chosenLeMethod == "gmres") { - LOG4CPLUS_INFO(logger, "Using GMRES method."); - // Transform the transition probability matrix to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); - if (precond == "ilu") { - gmm::gmres(*gmmA, x, b, gmm::ilu_precond>(*gmmA), 50, iter); - } else if (precond == "diagonal") { - gmm::gmres(*gmmA, x, b, gmm::diagonal_precond>(*gmmA), 50, iter); - } else if (precond == "ildlt") { - gmm::gmres(*gmmA, x, b, gmm::ildlt_precond>(*gmmA), 50, iter); - } else if (precond == "none") { - gmm::gmres(*gmmA, x, b, gmm::identity_matrix(), 50, iter); - } - - // Check if the solver converged and issue a warning otherwise. - if (iter.converged()) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); - } - delete gmmA; - } else if (chosenLeMethod == "jacobi") { - LOG4CPLUS_INFO(logger, "Using Jacobi method."); - uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b); - uint_fast64_t maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger(); - // Check if the solver converged and issue a warning otherwise. - if (iterations < maxIterations) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); - } - } - } + /*! + * Constructs a linear equation solver with the given parameters. + * + * @param method The method to use for linear equation solving. + * @param precision The precision to use for convergence detection. + * @param relative If set, the relative error rather than the absolute error is considered for convergence + * detection. + * @param maximalNumberOfIterations The maximal number of iterations do perform before iteration is aborted. + * @param preconditioner The preconditioner to use when solving linear equation systems. + * @param restart An optional argument that specifies after how many iterations restarted methods are + * supposed to actually to a restart. + */ + GmmxxLinearEquationSolver(SolutionMethod method, double precision, bool relative, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, uint_fast64_t restart = 0); - virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector* b, uint_fast64_t n = 1) const { - // Transform the transition probability A to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); - - // Set up some temporary variables so that we can just swap pointers instead of copying the result after each - // iteration. - std::vector* swap = nullptr; - std::vector* currentVector = &x; - std::vector* tmpVector = new std::vector(A.getRowCount()); - - // Now perform matrix-vector multiplication as long as we meet the bound. - for (uint_fast64_t i = 0; i < n; ++i) { - gmm::mult(*gmmxxMatrix, *currentVector, *tmpVector); - swap = tmpVector; - tmpVector = currentVector; - currentVector = swap; - - // If requested, add an offset to the current result vector. - if (b != nullptr) { - gmm::add(*b, *currentVector); - } - } - - // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, because - // the output is supposed to be stored in x. - if (n % 2 == 1) { - std::swap(x, *currentVector); - delete currentVector; - } else { - delete tmpVector; - } - - delete gmmxxMatrix; - } + /*! + * Constructs a linear equation solver with parameters being set according to the settings object. + */ + GmmxxLinearEquationSolver(); + + virtual AbstractLinearEquationSolver* clone() const override; + + virtual void solveEquationSystem(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; + + virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& A, std::vector& x, std::vector* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; private: /*! * Solves the linear equation system A*x = b given by the parameters using the Jacobi method. * * @param A The matrix specifying the coefficients of the linear equations. - * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but - * may be ignored. + * @param x The solution vector x. The initial values of x represent a guess of the real values to the + * solver, but may be set to zero. * @param b The right-hand side of the equation system. - * @returns The solution vector x of the system of linear equations as the content of the parameter x. - * @returns The number of iterations needed until convergence. + * @return The number of iterations needed until convergence if the solver converged and + * maximalNumberOfIteration otherwise. + * @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. + */ + uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const; + + /*! + * Retrieves the string representation of the solution method associated with this solver. + * + * @return The string representation of the solution method associated with this solver. + */ + std::string methodToString() const; + + /*! + * Retrieves the string representation of the preconditioner associated with this solver. + * + * @return The string representation of the preconditioner associated with this solver. */ - uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { - // Get the settings object to customize linear solving. - storm::settings::Settings* s = storm::settings::Settings::getInstance(); - - double precision = s->getOptionByLongName("precision").getArgument(0).getValueAsDouble(); - uint_fast64_t maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger(); - bool relative = s->getOptionByLongName("relative").getArgument(0).getValueAsBoolean(); - - // Get a Jacobi decomposition of the matrix A. - std::pair, storm::storage::SparseMatrix> jacobiDecomposition = A.getJacobiDecomposition(); - - // Convert the (inverted) diagonal matrix to gmm++'s format. - gmm::csr_matrix* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.second)); - // Convert the LU matrix to gmm++'s format. - gmm::csr_matrix* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.first)); - - LOG4CPLUS_INFO(logger, "Starting iterative Jacobi Solver."); - - // x_(k + 1) = D^-1 * (b - R * x_k) - // In x we keep a copy of the result for swapping in the loop (e.g. less copy-back). - std::vector* xNext = new std::vector(x.size()); - const std::vector* xCopy = xNext; - std::vector* xCurrent = &x; - - // Target vector for precision calculation. - std::vector tmpX(x.size()); - - // Set up additional environment variables. - uint_fast64_t iterationCount = 0; - bool converged = false; - - while (!converged && iterationCount < maxIterations) { - // R * x_k (xCurrent is x_k) -> xNext - gmm::mult(*gmmxxLU, *xCurrent, tmpX); - // b - R * x_k (xNext contains R * x_k) -> xNext - gmm::add(b, gmm::scaled(tmpX, -1.0), tmpX); - // D^-1 * (b - R * x_k) -> xNext - gmm::mult(*gmmxxDiagonalInverted, tmpX, *xNext); - - // Swap xNext with xCurrent so that the next iteration can use xCurrent again without having to copy the - // vector. - std::vector *const swap = xNext; - xNext = xCurrent; - xCurrent = swap; - - // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision(*xCurrent, *xNext, precision, relative); - - // Increase iteration count so we can abort if convergence is too slow. - ++iterationCount; - } - - // If the last iteration did not write to the original x we have to swap the contents, because the output has to - // be written to the parameter x. - if (xCurrent == xCopy) { - std::swap(x, *xCurrent); - } - - // As xCopy always points to the copy of x used for swapping, we can safely delete it. - delete xCopy; - - // Also delete the other dynamically allocated variables. - delete gmmxxDiagonalInverted; - delete gmmxxLU; - - return iterationCount; - } + std::string preconditionerToString() const; + // The method to use for solving linear equation systems. + SolutionMethod method; + + // The required precision for the iterative methods. + double precision; + + // Sets whether the relative or absolute error is to be considered for convergence detection. Not that this + // only applies to the Jacobi method for this solver. + bool relative; + + // The maximal number of iterations to do before iteration is aborted. + uint_fast64_t maximalNumberOfIterations; + + // The preconditioner to use when solving the linear equation system. + Preconditioner preconditioner; + + // A restart value that determines when restarted methods shall do so. + uint_fast64_t restart; }; - + } // namespace solver } // namespace storm diff --git a/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp b/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp new file mode 100644 index 000000000..74d976171 --- /dev/null +++ b/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp @@ -0,0 +1,11 @@ +#include "src/settings/Settings.h" + +bool GmmxxNondeterministicLinearEquationSolverOptionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* instance) -> bool { + instance->addOption(storm::settings::OptionBuilder("GmmxxNondeterminsticLinearEquationSolver", "maxiter", "i", "The maximal number of iterations to perform before iterative solving is aborted.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(10000).build()).build()); + + instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "precision", "", "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-6).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build()); + + instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "absolute", "", "Whether the relative or the absolute error is considered for deciding convergence.").build()); + + return true; +}); \ No newline at end of file diff --git a/src/solver/GmmxxNondeterministicLinearEquationSolver.h b/src/solver/GmmxxNondeterministicLinearEquationSolver.h index 205e724f1..33c5c97b6 100644 --- a/src/solver/GmmxxNondeterministicLinearEquationSolver.h +++ b/src/solver/GmmxxNondeterministicLinearEquationSolver.h @@ -22,7 +22,7 @@ namespace storm { virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& nondeterministicChoiceIndices, std::vector* b = nullptr, uint_fast64_t n = 1) const override { // Transform the transition probability matrix to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + std::unique_ptr> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); // Create vector for result of multiplication, which is reduced to the result vector after // each multiplication. @@ -42,14 +42,11 @@ namespace storm { storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices); } } - - // Delete intermediate results and return result. - delete gmmxxMatrix; } virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector const& nondeterministicChoiceIndices, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const override { // Transform the transition probability matrix to the gmm++ format to use its arithmetic. - gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + std::unique_ptr> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); // Set up the environment for the power method. bool multiplyResultMemoryProvided = true; @@ -102,7 +99,6 @@ namespace storm { } else if (!xMemoryProvided) { delete newX; } - delete gmmxxMatrix; if (!multiplyResultMemoryProvided) { delete multiplyResult;