From 023325b53d3edb3cddf20253af4c355ca7b930a0 Mon Sep 17 00:00:00 2001 From: dehnert Date: Fri, 17 Jun 2016 16:50:59 +0200 Subject: [PATCH] added tests for Eigen solver Former-commit-id: ede9efcee200c93b423cc89809c1027a447aa69c --- src/adapters/EigenAdapter.cpp | 1 + src/adapters/EigenAdapter.h | 36 ++++ src/adapters/GmmxxAdapter.h | 97 +++++----- src/settings/SettingsManager.cpp | 2 + .../modules/EigenEquationSolverSettings.cpp | 93 ++++++++++ .../modules/EigenEquationSolverSettings.h | 100 +++++++++++ src/settings/modules/MarkovChainSettings.cpp | 6 +- src/solver/EigenLinearEquationSolver.cpp | 117 ++++++++++++ src/solver/EigenLinearEquationSolver.h | 53 ++++++ src/solver/GmmxxLinearEquationSolver.h | 2 +- src/solver/SolverSelectionOptions.h | 2 +- src/utility/ConversionHelper.cpp | 5 - src/utility/ConversionHelper.h | 38 ---- src/utility/eigen.h | 3 + src/utility/solver.cpp | 7 + src/utility/solver.h | 6 + .../solver/EigenLinearEquationSolverTest.cpp | 170 ++++++++++++++++++ .../solver/GmmxxLinearEquationSolverTest.cpp | 4 +- 18 files changed, 640 insertions(+), 102 deletions(-) create mode 100644 src/adapters/EigenAdapter.cpp create mode 100644 src/adapters/EigenAdapter.h create mode 100644 src/settings/modules/EigenEquationSolverSettings.cpp create mode 100644 src/settings/modules/EigenEquationSolverSettings.h create mode 100644 src/solver/EigenLinearEquationSolver.cpp create mode 100644 src/solver/EigenLinearEquationSolver.h delete mode 100644 src/utility/ConversionHelper.cpp delete mode 100644 src/utility/ConversionHelper.h create mode 100644 src/utility/eigen.h create mode 100644 test/functional/solver/EigenLinearEquationSolverTest.cpp diff --git a/src/adapters/EigenAdapter.cpp b/src/adapters/EigenAdapter.cpp new file mode 100644 index 000000000..35aa9188c --- /dev/null +++ b/src/adapters/EigenAdapter.cpp @@ -0,0 +1 @@ +#include "src/adapters/EigenAdapter.h" \ No newline at end of file diff --git a/src/adapters/EigenAdapter.h b/src/adapters/EigenAdapter.h new file mode 100644 index 000000000..fc847094c --- /dev/null +++ b/src/adapters/EigenAdapter.h @@ -0,0 +1,36 @@ +#pragma once + +#include + +#include "src/utility/eigen.h" + +#include "src/storage/SparseMatrix.h" + +namespace storm { + namespace adapters { + + class EigenAdapter { + public: + /*! + * Converts a sparse matrix into a sparse matrix in the gmm++ format. + * @return A pointer to a row-major sparse matrix in gmm++ format. + */ + template + static std::unique_ptr> toEigenSparseMatrix(storm::storage::SparseMatrix const& matrix) { + std::vector> triplets; + triplets.reserve(matrix.getNonzeroEntryCount()); + + for (uint64_t row = 0; row < matrix.getRowCount(); ++row) { + for (auto const& element : matrix.getRow(row)) { + triplets.emplace_back(row, element.getColumn(), element.getValue()); + } + } + + std::unique_ptr> result = std::make_unique>(matrix.getRowCount(), matrix.getColumnCount()); + result->setFromTriplets(triplets.begin(), triplets.end()); + return result; + } + }; + + } +} \ No newline at end of file diff --git a/src/adapters/GmmxxAdapter.h b/src/adapters/GmmxxAdapter.h index ad4c82a0e..929a8303b 100644 --- a/src/adapters/GmmxxAdapter.h +++ b/src/adapters/GmmxxAdapter.h @@ -1,68 +1,59 @@ -/* - * GmmxxAdapter.h - * - * Created on: 24.12.2012 - * Author: Christian Dehnert - */ - #ifndef STORM_ADAPTERS_GMMXXADAPTER_H_ #define STORM_ADAPTERS_GMMXXADAPTER_H_ -#include -#include +#include +#include #include "src/utility/gmm.h" #include "src/storage/SparseMatrix.h" -#include "src/utility/ConversionHelper.h" #include "src/utility/macros.h" namespace storm { - -namespace adapters { - -class GmmxxAdapter { -public: - /*! - * Converts a sparse matrix into a sparse matrix in the gmm++ format. - * @return A pointer to a row-major sparse matrix in gmm++ format. - */ - template - static std::unique_ptr> toGmmxxSparseMatrix(storm::storage::SparseMatrix const& matrix) { - uint_fast64_t realNonZeros = matrix.getEntryCount(); - STORM_LOG_DEBUG("Converting matrix with " << realNonZeros << " non-zeros to gmm++ format."); - - // Prepare the resulting matrix. - 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()); - - // Copy columns and values. - std::vector values; - values.reserve(matrix.getEntryCount()); + + namespace adapters { - // To match the correct vector type for gmm, we create the vector with the exact same type. - decltype(result->ir) columns; - columns.reserve(matrix.getEntryCount()); + class GmmxxAdapter { + public: + /*! + * Converts a sparse matrix into a sparse matrix in the gmm++ format. + * @return A pointer to a row-major sparse matrix in gmm++ format. + */ + template + static std::unique_ptr> toGmmxxSparseMatrix(storm::storage::SparseMatrix const& matrix) { + uint_fast64_t realNonZeros = matrix.getEntryCount(); + STORM_LOG_DEBUG("Converting matrix with " << realNonZeros << " non-zeros to gmm++ format."); + + // Prepare the resulting matrix. + 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()); + + // Copy columns and values. + std::vector values; + values.reserve(matrix.getEntryCount()); + + // To match the correct vector type for gmm, we create the vector with the exact same type. + decltype(result->ir) columns; + columns.reserve(matrix.getEntryCount()); + + for (auto const& entry : matrix) { + columns.emplace_back(entry.getColumn()); + values.emplace_back(entry.getValue()); + } + + std::swap(result->ir, columns); + std::swap(result->pr, values); + + STORM_LOG_DEBUG("Done converting matrix to gmm++ format."); + + return result; + } + }; - for (auto const& entry : matrix) { - columns.emplace_back(entry.getColumn()); - values.emplace_back(entry.getValue()); - } - - std::swap(result->ir, columns); - std::swap(result->pr, values); - - STORM_LOG_DEBUG("Done converting matrix to gmm++ format."); - - return result; - } -}; - -} //namespace adapters - -} //namespace storm + } // namespace adapters +} // namespace storm #endif /* STORM_ADAPTERS_GMMXXADAPTER_H_ */ diff --git a/src/settings/SettingsManager.cpp b/src/settings/SettingsManager.cpp index 8e79c5457..452b20fd3 100644 --- a/src/settings/SettingsManager.cpp +++ b/src/settings/SettingsManager.cpp @@ -20,6 +20,7 @@ #include "src/settings/modules/CounterexampleGeneratorSettings.h" #include "src/settings/modules/CuddSettings.h" #include "src/settings/modules/SylvanSettings.h" +#include "src/settings/modules/EigenEquationSolverSettings.h" #include "src/settings/modules/GmmxxEquationSolverSettings.h" #include "src/settings/modules/NativeEquationSolverSettings.h" #include "src/settings/modules/BisimulationSettings.h" @@ -508,6 +509,7 @@ namespace storm { storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); + storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); diff --git a/src/settings/modules/EigenEquationSolverSettings.cpp b/src/settings/modules/EigenEquationSolverSettings.cpp new file mode 100644 index 000000000..568650c0a --- /dev/null +++ b/src/settings/modules/EigenEquationSolverSettings.cpp @@ -0,0 +1,93 @@ +#include "src/settings/modules/EigenEquationSolverSettings.h" + +#include "src/settings/Option.h" +#include "src/settings/OptionBuilder.h" +#include "src/settings/ArgumentBuilder.h" +#include "src/settings/Argument.h" + +#include "src/settings/SettingsManager.h" +#include "src/settings/modules/MarkovChainSettings.h" +#include "src/solver/SolverSelectionOptions.h" + +namespace storm { + namespace settings { + namespace modules { + + const std::string EigenEquationSolverSettings::moduleName = "eigen"; + const std::string EigenEquationSolverSettings::techniqueOptionName = "method"; + const std::string EigenEquationSolverSettings::preconditionOptionName = "precond"; + const std::string EigenEquationSolverSettings::maximalIterationsOptionName = "maxiter"; + const std::string EigenEquationSolverSettings::maximalIterationsOptionShortName = "i"; + const std::string EigenEquationSolverSettings::precisionOptionName = "precision"; + + EigenEquationSolverSettings::EigenEquationSolverSettings() : ModuleSettings(moduleName) { + std::vector methods = {"sparselu", "bicgstab"}; + this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the eigen solver. Available are {sparselu, bicgstab}.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("sparselu").build()).build()); + + // Register available preconditioners. + std::vector preconditioner = {"ilu", "diagonal", "none"}; + this->addOption(storm::settings::OptionBuilder(moduleName, preconditionOptionName, true, "The preconditioning technique used for solving linear equation systems. 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()); + + this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build()); + + this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "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()); + } + + bool EigenEquationSolverSettings::isLinearEquationSystemMethodSet() const { + return this->getOption(techniqueOptionName).getHasOptionBeenSet(); + } + + EigenEquationSolverSettings::LinearEquationMethod EigenEquationSolverSettings::getLinearEquationSystemMethod() const { + std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString(); + if (linearEquationSystemTechniqueAsString == "sparselu") { + return EigenEquationSolverSettings::LinearEquationMethod::SparseLU; + } else if (linearEquationSystemTechniqueAsString == "bicgstab") { + return EigenEquationSolverSettings::LinearEquationMethod::Bicgstab; + } + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected."); + } + + bool EigenEquationSolverSettings::isPreconditioningMethodSet() const { + return this->getOption(preconditionOptionName).getHasOptionBeenSet(); + } + + EigenEquationSolverSettings::PreconditioningMethod EigenEquationSolverSettings::getPreconditioningMethod() const { + std::string PreconditioningMethodAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString(); + if (PreconditioningMethodAsString == "ilu") { + return EigenEquationSolverSettings::PreconditioningMethod::Ilu; + } else if (PreconditioningMethodAsString == "diagonal") { + return EigenEquationSolverSettings::PreconditioningMethod::Diagonal; + } else if (PreconditioningMethodAsString == "none") { + return EigenEquationSolverSettings::PreconditioningMethod::None; + } + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << PreconditioningMethodAsString << "' selected."); + } + + bool EigenEquationSolverSettings::isMaximalIterationCountSet() const { + return this->getOption(maximalIterationsOptionName).getHasOptionBeenSet(); + } + + uint_fast64_t EigenEquationSolverSettings::getMaximalIterationCount() const { + return this->getOption(maximalIterationsOptionName).getArgumentByName("count").getValueAsUnsignedInteger(); + } + + bool EigenEquationSolverSettings::isPrecisionSet() const { + return this->getOption(precisionOptionName).getHasOptionBeenSet(); + } + + double EigenEquationSolverSettings::getPrecision() const { + return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble(); + } + + bool EigenEquationSolverSettings::check() const { + // This list does not include the precision, because this option is shared with other modules. + bool optionsSet = isLinearEquationSystemMethodSet() || isPreconditioningMethodSet() || isMaximalIterationCountSet(); + + STORM_LOG_WARN_COND(storm::settings::getModule().getEquationSolver() == storm::solver::EquationSolverType::Gmmxx || !optionsSet, "eigen is not selected as the preferred equation solver, so setting options for eigen might have no effect."); + + return true; + } + + } // namespace modules + } // namespace settings +} // namespace storm diff --git a/src/settings/modules/EigenEquationSolverSettings.h b/src/settings/modules/EigenEquationSolverSettings.h new file mode 100644 index 000000000..d415b2c31 --- /dev/null +++ b/src/settings/modules/EigenEquationSolverSettings.h @@ -0,0 +1,100 @@ +#pragma once + +#include "src/settings/modules/ModuleSettings.h" + +namespace storm { + namespace settings { + namespace modules { + + /*! + * This class represents the settings for Eigen. + */ + class EigenEquationSolverSettings : public ModuleSettings { + public: + // An enumeration of all available methods for solving linear equations. + enum class LinearEquationMethod { SparseLU, Bicgstab }; + + // An enumeration of all available preconditioning methods. + enum class PreconditioningMethod { Ilu, Diagonal, None }; + + // An enumeration of all available convergence criteria. + enum class ConvergenceCriterion { Absolute, Relative }; + + /*! + * Creates a new set of Eigen settings. + */ + EigenEquationSolverSettings(); + + /*! + * Retrieves whether the linear equation system method has been set. + * + * @return True iff the linear equation system method has been set. + */ + bool isLinearEquationSystemMethodSet() const; + + /*! + * Retrieves the method that is to be used for solving systems of linear equations. + * + * @return The method to use. + */ + LinearEquationMethod getLinearEquationSystemMethod() const; + + /*! + * Retrieves whether the preconditioning method has been set. + * + * @return True iff the preconditioning method has been set. + */ + bool isPreconditioningMethodSet() const; + + /*! + * Retrieves the method that is to be used for preconditioning solving systems of linear equations. + * + * @return The method to use. + */ + PreconditioningMethod getPreconditioningMethod() const; + + /*! + * Retrieves whether the maximal iteration count has been set. + * + * @return True iff the maximal iteration count has been set. + */ + bool isMaximalIterationCountSet() const; + + /*! + * Retrieves the maximal number of iterations to perform until giving up on converging. + * + * @return The maximal iteration count. + */ + uint_fast64_t getMaximalIterationCount() const; + + /*! + * Retrieves whether the precision has been set. + * + * @return True iff the precision has been set. + */ + bool isPrecisionSet() const; + + /*! + * Retrieves the precision that is used for detecting convergence. + * + * @return The precision to use for detecting convergence. + */ + double getPrecision() const; + + bool check() const override; + + // The name of the module. + static const std::string moduleName; + + private: + // Define the string names of the options as constants. + static const std::string techniqueOptionName; + static const std::string preconditionOptionName; + static const std::string maximalIterationsOptionName; + static const std::string maximalIterationsOptionShortName; + static const std::string precisionOptionName; + }; + + } // namespace modules + } // namespace settings +} // namespace storm diff --git a/src/settings/modules/MarkovChainSettings.cpp b/src/settings/modules/MarkovChainSettings.cpp index d436463a4..589637103 100644 --- a/src/settings/modules/MarkovChainSettings.cpp +++ b/src/settings/modules/MarkovChainSettings.cpp @@ -41,9 +41,9 @@ namespace storm { this->addOption(storm::settings::OptionBuilder(moduleName, engineOptionName, false, "Sets which engine is used for model building and model checking.").setShortName(engineOptionShortName) .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the engine to use. Available are {sparse, hybrid, dd, expl, abs}.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(engines)).setDefaultValueString("sparse").build()).build()); - std::vector linearEquationSolver = {"gmm++", "native"}; + std::vector linearEquationSolver = {"gmm++", "native", "eigen"}; this->addOption(storm::settings::OptionBuilder(moduleName, eqSolverOptionName, false, "Sets which solver is preferred for solving systems of linear equations.") - .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the solver to prefer. Available are: gmm++ and native.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(linearEquationSolver)).setDefaultValueString("gmm++").build()).build()); + .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the solver to prefer. Available are: gmm++, native and eigen.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(linearEquationSolver)).setDefaultValueString("gmm++").build()).build()); std::vector ddLibraries = {"cudd", "sylvan"}; this->addOption(storm::settings::OptionBuilder(moduleName, ddLibraryOptionName, false, "Sets which library is preferred for decision-diagram operations.") @@ -85,6 +85,8 @@ namespace storm { return storm::solver::EquationSolverType::Gmmxx; } else if (equationSolverName == "native") { return storm::solver::EquationSolverType::Native; + } else if (equationSolverName == "eigen") { + return storm::solver::EquationSolverType::Eigen; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown equation solver '" << equationSolverName << "'."); } diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp new file mode 100644 index 000000000..7a6c83231 --- /dev/null +++ b/src/solver/EigenLinearEquationSolver.cpp @@ -0,0 +1,117 @@ +#include "src/solver/EigenLinearEquationSolver.h" + +#include "src/adapters/EigenAdapter.h" + +#include "src/settings/SettingsManager.h" +#include "src/settings/modules/EigenEquationSolverSettings.h" + +namespace storm { + namespace solver { + + template + EigenLinearEquationSolver::EigenLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix(A)), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner) { + // Intentionally left empty. + } + + template + EigenLinearEquationSolver::EigenLinearEquationSolver(storm::storage::SparseMatrix const& A) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix(A)) { + // Get the settings object to customize linear solving. + storm::settings::modules::EigenEquationSolverSettings const& settings = storm::settings::getModule(); + + // Get appropriate settings. + maximalNumberOfIterations = settings.getMaximalIterationCount(); + precision = settings.getPrecision(); + + // Determine the method to be used. + storm::settings::modules::EigenEquationSolverSettings::LinearEquationMethod methodAsSetting = settings.getLinearEquationSystemMethod(); + if (methodAsSetting == storm::settings::modules::EigenEquationSolverSettings::LinearEquationMethod::Bicgstab) { + method = SolutionMethod::Bicgstab; + } else if (methodAsSetting == storm::settings::modules::EigenEquationSolverSettings::LinearEquationMethod::SparseLU) { + method = SolutionMethod::SparseLU; + } + + // Check which preconditioner to use. + storm::settings::modules::EigenEquationSolverSettings::PreconditioningMethod preconditionAsSetting = settings.getPreconditioningMethod(); + if (preconditionAsSetting == storm::settings::modules::EigenEquationSolverSettings::PreconditioningMethod::Ilu) { + preconditioner = Preconditioner::Ilu; + } else if (preconditionAsSetting == storm::settings::modules::EigenEquationSolverSettings::PreconditioningMethod::Diagonal) { + preconditioner = Preconditioner::Diagonal; + } else if (preconditionAsSetting == storm::settings::modules::EigenEquationSolverSettings::PreconditioningMethod::None) { + preconditioner = Preconditioner::None; + } + } + + template + void EigenLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + // Translate the vectors x and b into Eigen's format. + Eigen::VectorXd eigenB(b.size()); + for (uint64_t index = 0; index < b.size(); ++index) { + eigenB[index] = b[index]; + } + + Eigen::VectorXd eigenX(x.size()); + for (uint64_t index = 0; index < x.size(); ++index) { + eigenX[index] = x[index]; + } + + if (method == SolutionMethod::SparseLU) { + Eigen::SparseLU, Eigen::COLAMDOrdering> solver; + solver.compute(*eigenA); + solver._solve(eigenB, eigenX); + } else { + if (preconditioner == Preconditioner::Ilu) { + Eigen::BiCGSTAB, Eigen::IncompleteLUT> solver; + solver.compute(*eigenA); + solver._solve(eigenB, eigenX); + } else if (preconditioner == Preconditioner::Diagonal) { + Eigen::BiCGSTAB, Eigen::DiagonalPreconditioner> solver; + solver.compute(*eigenA); + solver._solve(eigenB, eigenX); + } else { + Eigen::BiCGSTAB, Eigen::IdentityPreconditioner> solver; + solver.compute(*eigenA); + solver._solve(eigenB, eigenX); + } + } + + // Translate the solution from Eigen's format into our representation. + for (uint64_t index = 0; index < eigenX.size(); ++index) { + x[index] = eigenX[index]; + } + } + + template + void EigenLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { + + // Translate the vectors x and b into Eigen's format. + std::unique_ptr eigenB; + if (b != nullptr) { + eigenB = std::make_unique(b->size()); + for (uint64_t index = 0; index < b->size(); ++index) { + (*eigenB)[index] = (*b)[index]; + } + } + Eigen::VectorXd eigenX(x.size()); + for (uint64_t index = 0; index < x.size(); ++index) { + eigenX[index] = x[index]; + } + + // Perform n matrix-vector multiplications. + for (uint64_t iteration = 0; iteration < n; ++iteration) { + eigenX = *eigenA * eigenX; + if (eigenB != nullptr) { + eigenX += *eigenB; + } + } + + // Translate the solution from Eigen's format into our representation. + for (uint64_t index = 0; index < eigenX.size(); ++index) { + x[index] = eigenX[index]; + } + } + + template class EigenLinearEquationSolver; + + + } +} \ No newline at end of file diff --git a/src/solver/EigenLinearEquationSolver.h b/src/solver/EigenLinearEquationSolver.h new file mode 100644 index 000000000..91e84f8f0 --- /dev/null +++ b/src/solver/EigenLinearEquationSolver.h @@ -0,0 +1,53 @@ +#pragma once + +#include "src/solver/LinearEquationSolver.h" + +#include "src/utility/eigen.h" + +namespace storm { + namespace solver { + + /*! + * A class that uses the Eigen library to implement the LinearEquationSolver interface. + */ + template + class EigenLinearEquationSolver : public LinearEquationSolver { + public: + enum class SolutionMethod { + SparseLU, Bicgstab + }; + + enum class Preconditioner { + Ilu, Diagonal, None + }; + + EigenLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner); + + EigenLinearEquationSolver(storm::storage::SparseMatrix const& A); + + virtual void 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; + + private: + // A pointer to the original sparse matrix given to this solver. + storm::storage::SparseMatrix const* originalA; + + // The (eigen) matrix associated with this equation solver. + std::unique_ptr> eigenA; + + // The method to use for solving linear equation systems. + SolutionMethod method; + + // The required precision for the iterative methods. + double precision; + + // 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; + }; + + } +} \ No newline at end of file diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index a8a029728..a6c402f89 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -97,7 +97,7 @@ namespace storm { // The preconditioner to use when solving the linear equation system. Preconditioner preconditioner; - // Sets whether the relative or absolute error is to be considered for convergence detection. Not that this + // Sets whether the relative or absolute error is to be considered for convergence detection. Note that this // only applies to the Jacobi method for this solver. bool relative; diff --git a/src/solver/SolverSelectionOptions.h b/src/solver/SolverSelectionOptions.h index 1379343c9..49f28d1de 100644 --- a/src/solver/SolverSelectionOptions.h +++ b/src/solver/SolverSelectionOptions.h @@ -9,7 +9,7 @@ namespace storm { ExtendEnumsWithSelectionField(MinMaxTechnique, PolicyIteration, ValueIteration) ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk) - ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Topological) + ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Topological) ExtendEnumsWithSelectionField(SmtSolverType, Z3, Mathsat) } } diff --git a/src/utility/ConversionHelper.cpp b/src/utility/ConversionHelper.cpp deleted file mode 100644 index cb87b44e1..000000000 --- a/src/utility/ConversionHelper.cpp +++ /dev/null @@ -1,5 +0,0 @@ -#include "ConversionHelper.h" - -std::vector>* storm::utility::ConversionHelper::toUnsignedLongLong(std::vector>* vectorPtr) { - return reinterpret_cast> *>(vectorPtr); -} \ No newline at end of file diff --git a/src/utility/ConversionHelper.h b/src/utility/ConversionHelper.h deleted file mode 100644 index b87b520aa..000000000 --- a/src/utility/ConversionHelper.h +++ /dev/null @@ -1,38 +0,0 @@ -/* - * ConversionHelper.h - * - * Created on: 14.09.2013 - * Author: Philipp Berger - * - * WARNING: This file REQUIRES -no-strict-aliasing! - */ - -#ifndef STORM_UTILITY_CONVERSIONHELPER_H_ -#define STORM_UTILITY_CONVERSIONHELPER_H_ - -#include -#include -#include - -static_assert(sizeof(unsigned long long) == sizeof(uint_fast64_t), "This program uses the GMM Backend and therefor requires unsigned long long and uint_fast64_t to be of the same size!"); - - -namespace storm { - namespace utility { - - class ConversionHelper { - public: - /*! - * Converts a pointer to a std::vector to std::vector - */ - static std::vector>* toUnsignedLongLong(std::vector>* vectorPtr); - - private: - ConversionHelper() {} - ConversionHelper(ConversionHelper& other) {} - ~ConversionHelper() {} - }; - } -} - -#endif // STORM_UTILITY_CONVERSIONHELPER_H_ \ No newline at end of file diff --git a/src/utility/eigen.h b/src/utility/eigen.h new file mode 100644 index 000000000..d3816b788 --- /dev/null +++ b/src/utility/eigen.h @@ -0,0 +1,3 @@ +#pragma once + +#include diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp index 060fef4c9..6c1d56077 100644 --- a/src/utility/solver.cpp +++ b/src/utility/solver.cpp @@ -7,6 +7,7 @@ #include "src/solver/SymbolicGameSolver.h" #include "src/solver/NativeLinearEquationSolver.h" #include "src/solver/GmmxxLinearEquationSolver.h" +#include "src/solver/EigenLinearEquationSolver.h" #include "src/solver/GameSolver.h" #include "src/solver/NativeMinMaxLinearEquationSolver.h" @@ -50,6 +51,7 @@ namespace storm { switch (equationSolver) { case storm::solver::EquationSolverType::Gmmxx: return std::unique_ptr>(new storm::solver::GmmxxLinearEquationSolver(matrix)); case storm::solver::EquationSolverType::Native: return std::unique_ptr>(new storm::solver::NativeLinearEquationSolver(matrix)); + case storm::solver::EquationSolverType::Eigen: return std::unique_ptr>(new storm::solver::EigenLinearEquationSolver(matrix)); default: return std::unique_ptr>(new storm::solver::GmmxxLinearEquationSolver(matrix)); } } @@ -58,6 +60,11 @@ namespace storm { std::unique_ptr> GmmxxLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { return std::unique_ptr>(new storm::solver::GmmxxLinearEquationSolver(matrix)); } + + template + std::unique_ptr> EigenLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { + return std::unique_ptr>(new storm::solver::EigenLinearEquationSolver(matrix)); + } template NativeLinearEquationSolverFactory::NativeLinearEquationSolverFactory() { diff --git a/src/utility/solver.h b/src/utility/solver.h index ef8120b65..c0770605f 100644 --- a/src/utility/solver.h +++ b/src/utility/solver.h @@ -106,6 +106,12 @@ namespace storm { virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const override; }; + template + class EigenLinearEquationSolverFactory : public LinearEquationSolverFactory { + public: + virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const override; + }; + template class MinMaxLinearEquationSolverFactory { public: diff --git a/test/functional/solver/EigenLinearEquationSolverTest.cpp b/test/functional/solver/EigenLinearEquationSolverTest.cpp new file mode 100644 index 000000000..87703de9a --- /dev/null +++ b/test/functional/solver/EigenLinearEquationSolverTest.cpp @@ -0,0 +1,170 @@ +#include "gtest/gtest.h" +#include "storm-config.h" + +#include "src/solver/EigenLinearEquationSolver.h" +#include "src/settings/SettingsManager.h" + +#include "src/settings/modules/GmmxxEquationSolverSettings.h" + +TEST(EigenLinearEquationSolver, SolveWithStandardOptions) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder builder; + ASSERT_NO_THROW(builder.addNextValue(0, 0, 2)); + ASSERT_NO_THROW(builder.addNextValue(0, 1, 4)); + ASSERT_NO_THROW(builder.addNextValue(0, 2, -2)); + ASSERT_NO_THROW(builder.addNextValue(1, 0, 4)); + ASSERT_NO_THROW(builder.addNextValue(1, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(1, 2, 5)); + ASSERT_NO_THROW(builder.addNextValue(2, 0, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 2, 3)); + + storm::storage::SparseMatrix A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(3); + std::vector b = {16, -4, -7}; + + ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver solver(A)); + + storm::solver::EigenLinearEquationSolver solver(A); + ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); +} + +TEST(EigenLinearEquationSolver, SparseLU) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder builder; + ASSERT_NO_THROW(builder.addNextValue(0, 0, 2)); + ASSERT_NO_THROW(builder.addNextValue(0, 1, 4)); + ASSERT_NO_THROW(builder.addNextValue(0, 2, -2)); + ASSERT_NO_THROW(builder.addNextValue(1, 0, 4)); + ASSERT_NO_THROW(builder.addNextValue(1, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(1, 2, 5)); + ASSERT_NO_THROW(builder.addNextValue(2, 0, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 2, 3)); + + storm::storage::SparseMatrix A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(3); + std::vector b = {16, -4, -7}; + + ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::SparseLU, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::None)); + + storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::SparseLU, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::None); + ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); +} + +TEST(EigenLinearEquationSolver, BiCGSTAB) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder builder; + ASSERT_NO_THROW(builder.addNextValue(0, 0, 2)); + ASSERT_NO_THROW(builder.addNextValue(0, 1, 4)); + ASSERT_NO_THROW(builder.addNextValue(0, 2, -2)); + ASSERT_NO_THROW(builder.addNextValue(1, 0, 4)); + ASSERT_NO_THROW(builder.addNextValue(1, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(1, 2, 5)); + ASSERT_NO_THROW(builder.addNextValue(2, 0, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 2, 3)); + + storm::storage::SparseMatrix A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(3); + std::vector b = {16, -4, -7}; + + ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::None)); + + storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::None); + ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); +} + +TEST(EigenLinearEquationSolver, BiCGSTAB_Ilu) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder builder; + ASSERT_NO_THROW(builder.addNextValue(0, 0, 2)); + ASSERT_NO_THROW(builder.addNextValue(0, 1, 4)); + ASSERT_NO_THROW(builder.addNextValue(0, 2, -2)); + ASSERT_NO_THROW(builder.addNextValue(1, 0, 4)); + ASSERT_NO_THROW(builder.addNextValue(1, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(1, 2, 5)); + ASSERT_NO_THROW(builder.addNextValue(2, 0, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 2, 3)); + + storm::storage::SparseMatrix A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(3); + std::vector b = {16, -4, -7}; + + ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::Ilu)); + + storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::Ilu); + ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); +} + +TEST(EigenLinearEquationSolver, BiCGSTAB_Diagonal) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder builder; + ASSERT_NO_THROW(builder.addNextValue(0, 0, 2)); + ASSERT_NO_THROW(builder.addNextValue(0, 1, 4)); + ASSERT_NO_THROW(builder.addNextValue(0, 2, -2)); + ASSERT_NO_THROW(builder.addNextValue(1, 0, 4)); + ASSERT_NO_THROW(builder.addNextValue(1, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(1, 2, 5)); + ASSERT_NO_THROW(builder.addNextValue(2, 0, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 1, -1)); + ASSERT_NO_THROW(builder.addNextValue(2, 2, 3)); + + storm::storage::SparseMatrix A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(3); + std::vector b = {16, -4, -7}; + + ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::Diagonal)); + + storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::Diagonal); + ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); +} + +TEST(EigenLinearEquationSolver, MatrixVectorMultplication) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder builder; + ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5)); + ASSERT_NO_THROW(builder.addNextValue(0, 4, 0.5)); + ASSERT_NO_THROW(builder.addNextValue(1, 2, 0.5)); + ASSERT_NO_THROW(builder.addNextValue(1, 4, 0.5)); + ASSERT_NO_THROW(builder.addNextValue(2, 3, 0.5)); + ASSERT_NO_THROW(builder.addNextValue(2, 4, 0.5)); + ASSERT_NO_THROW(builder.addNextValue(3, 4, 1)); + ASSERT_NO_THROW(builder.addNextValue(4, 4, 1)); + + storm::storage::SparseMatrix A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(5); + x[4] = 1; + + storm::solver::EigenLinearEquationSolver solver(A); + ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(x, nullptr, 4)); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); +} \ No newline at end of file diff --git a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp index 21e74fe7c..2820678d2 100644 --- a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp +++ b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp @@ -167,7 +167,7 @@ TEST(GmmxxLinearEquationSolver, gmresilu) { ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver solver(A, storm::solver::GmmxxLinearEquationSolver::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver::Preconditioner::Ilu, true, 50)); - storm::solver::GmmxxLinearEquationSolver solver(A, storm::solver::GmmxxLinearEquationSolver::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver::Preconditioner::None, true, 50); + storm::solver::GmmxxLinearEquationSolver solver(A, storm::solver::GmmxxLinearEquationSolver::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver::Preconditioner::Ilu, true, 50); ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); @@ -195,7 +195,7 @@ TEST(GmmxxLinearEquationSolver, gmresdiag) { ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver solver(A, storm::solver::GmmxxLinearEquationSolver::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver::Preconditioner::Diagonal, true, 50)); - storm::solver::GmmxxLinearEquationSolver solver(A, storm::solver::GmmxxLinearEquationSolver::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver::Preconditioner::None, true, 50); + storm::solver::GmmxxLinearEquationSolver solver(A, storm::solver::GmmxxLinearEquationSolver::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver::Preconditioner::Diagonal, true, 50); ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision());