From 48b82e7b14a91978eea99c6a1f534ef22fc69c71 Mon Sep 17 00:00:00 2001 From: TimQu Date: Tue, 15 Nov 2016 19:28:48 +0100 Subject: [PATCH] refactored auxData in linear equation solvers Former-commit-id: 9e0deb255caf6e6292409e0e85178ff72dd392fd --- .../csl/helper/SparseCtmcCslHelper.cpp | 1 - .../pcaa/SparseMaPcaaWeightVectorChecker.cpp | 1 - src/solver/GmmxxLinearEquationSolver.cpp | 141 +++++++----------- src/solver/GmmxxLinearEquationSolver.h | 24 ++- src/solver/LinearEquationSolver.cpp | 58 ++----- src/solver/LinearEquationSolver.h | 41 +---- src/solver/MinMaxLinearEquationSolver.h | 2 +- src/solver/NativeLinearEquationSolver.cpp | 106 ++++--------- src/solver/NativeLinearEquationSolver.h | 14 +- .../StandardMinMaxLinearEquationSolver.cpp | 50 ++++--- .../StandardMinMaxLinearEquationSolver.h | 10 +- src/utility/policyguessing.cpp | 4 +- .../solver/GmmxxLinearEquationSolverTest.cpp | 78 ++++++---- 13 files changed, 198 insertions(+), 332 deletions(-) diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp index c51563507..df8e06c4a 100644 --- a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp +++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp @@ -629,7 +629,6 @@ namespace storm { } std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(uniformizedMatrix)); - solver->allocateAuxMemory(storm::solver::LinearEquationSolverOperation::MultiplyRepeatedly); if (!useMixedPoissonProbabilities && std::get<0>(foxGlynnResult) > 1) { // Perform the matrix-vector multiplications (without adding). diff --git a/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp b/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp index 728736d34..2e4e13439 100644 --- a/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp +++ b/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp @@ -365,7 +365,6 @@ namespace storm { storm::storage::SparseMatrix linEqMatrix = PS.toPS.selectRowsFromRowGroups(optimalChoicesAtCurrentEpoch, true); linEqMatrix.convertToEquationSystem(); linEq.solver = linEq.factory.create(std::move(linEqMatrix)); - linEq.solver->allocateAuxMemory(storm::solver::LinearEquationSolverOperation::SolveEquations); } // Get the results for the individual objectives. diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index fd15e7e27..f69eb7fd3 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -111,12 +111,12 @@ namespace storm { } template - GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiMemory(nullptr) { + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings) { this->setMatrix(A); } template - GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiMemory(nullptr) { + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings) { this->setMatrix(std::move(A)); } @@ -124,66 +124,79 @@ namespace storm { void GmmxxLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { localA.reset(); this->A = &A; - gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); - jacobiDecomposition.reset(); + resetAuxiliaryData(); } template void GmmxxLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { localA = std::make_unique>(std::move(A)); this->A = localA.get(); - gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(*localA); + resetAuxiliaryData(); } template bool GmmxxLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { auto method = this->getSettings().getSolutionMethod(); auto preconditioner = this->getSettings().getPreconditioner(); - STORM_LOG_INFO("Using method '" << method << "' with preconditioner '" << preconditioner << "' (max. " << this->getSettings().getMaximalNumberOfIterations() << " iterations)."); + STORM_LOG_DEBUG("Using method '" << method << "' with preconditioner '" << preconditioner << "' (max. " << this->getSettings().getMaximalNumberOfIterations() << " iterations)."); if (method == GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi && preconditioner != GmmxxLinearEquationSolverSettings::Preconditioner::None) { STORM_LOG_WARN("Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored."); } if (method == GmmxxLinearEquationSolverSettings::SolutionMethod::Bicgstab || method == GmmxxLinearEquationSolverSettings::SolutionMethod::Qmr || method == GmmxxLinearEquationSolverSettings::SolutionMethod::Gmres) { + + // Translate the matrix into gmm++ format (if not already done) + if(!gmmxxA) { + gmmxxA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(*A); + } + + // Make sure that the requested preconditioner is available + if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::Ilu && !iluPreconditioner) { + iluPreconditioner = std::make_unique>>(*gmmxxA); + } else if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::Diagonal) { + diagonalPreconditioner = std::make_unique>>(*gmmxxA); + } + // Prepare an iteration object that determines the accuracy and the maximum number of iterations. gmm::iteration iter(this->getSettings().getPrecision(), 0, this->getSettings().getMaximalNumberOfIterations()); + // Invoke gmm with the corresponding settings if (method == GmmxxLinearEquationSolverSettings::SolutionMethod::Bicgstab) { if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::Ilu) { - gmm::bicgstab(*gmmxxMatrix, x, b, gmm::ilu_precond>(*gmmxxMatrix), iter); + gmm::bicgstab(*gmmxxA, x, b, *iluPreconditioner, iter); } else if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::Diagonal) { - gmm::bicgstab(*gmmxxMatrix, x, b, gmm::diagonal_precond>(*gmmxxMatrix), iter); + gmm::bicgstab(*gmmxxA, x, b, *diagonalPreconditioner, iter); } else if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::None) { - gmm::bicgstab(*gmmxxMatrix, x, b, gmm::identity_matrix(), iter); + gmm::bicgstab(*gmmxxA, x, b, gmm::identity_matrix(), iter); } } else if (method == GmmxxLinearEquationSolverSettings::SolutionMethod::Qmr) { if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::Ilu) { - gmm::qmr(*gmmxxMatrix, x, b, gmm::ilu_precond>(*gmmxxMatrix), iter); + gmm::qmr(*gmmxxA, x, b, *iluPreconditioner, iter); } else if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::Diagonal) { - gmm::qmr(*gmmxxMatrix, x, b, gmm::diagonal_precond>(*gmmxxMatrix), iter); + gmm::qmr(*gmmxxA, x, b, *diagonalPreconditioner, iter); } else if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::None) { - gmm::qmr(*gmmxxMatrix, x, b, gmm::identity_matrix(), iter); + gmm::qmr(*gmmxxA, x, b, gmm::identity_matrix(), iter); } } else if (method == GmmxxLinearEquationSolverSettings::SolutionMethod::Gmres) { if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::Ilu) { - gmm::gmres(*gmmxxMatrix, x, b, gmm::ilu_precond>(*gmmxxMatrix), this->getSettings().getNumberOfIterationsUntilRestart(), iter); + gmm::gmres(*gmmxxA, x, b, *iluPreconditioner, this->getSettings().getNumberOfIterationsUntilRestart(), iter); } else if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::Diagonal) { - gmm::gmres(*gmmxxMatrix, x, b, gmm::diagonal_precond>(*gmmxxMatrix), this->getSettings().getNumberOfIterationsUntilRestart(), iter); + gmm::gmres(*gmmxxA, x, b, *diagonalPreconditioner, this->getSettings().getNumberOfIterationsUntilRestart(), iter); } else if (preconditioner == GmmxxLinearEquationSolverSettings::Preconditioner::None) { - gmm::gmres(*gmmxxMatrix, x, b, gmm::identity_matrix(), this->getSettings().getNumberOfIterationsUntilRestart(), iter); + gmm::gmres(*gmmxxA, x, b, gmm::identity_matrix(), this->getSettings().getNumberOfIterationsUntilRestart(), iter); } } // Check if the solver converged and issue a warning otherwise. if (iter.converged()) { - STORM_LOG_INFO("Iterative solver converged after " << iter.get_iteration() << " iterations."); + STORM_LOG_DEBUG("Iterative solver converged after " << iter.get_iteration() << " iterations."); return true; } else { STORM_LOG_WARN("Iterative solver did not converge."); return false; } } else if (method == GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi) { - uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(*A, x, b); + uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(x, b); // Check if the solver converged and issue a warning otherwise. if (iterations < this->getSettings().getMaximalNumberOfIterations()) { @@ -200,39 +213,44 @@ namespace storm { template void GmmxxLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { + if(!gmmxxA) { + gmmxxA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(*A); + } if (b) { - gmm::mult_add(*gmmxxMatrix, x, *b, result); + gmm::mult_add(*gmmxxA, x, *b, result); } else { - gmm::mult(*gmmxxMatrix, x, result); + gmm::mult(*gmmxxA, x, result); } } template - uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { - bool allocatedAuxMemory = !this->hasAuxMemory(LinearEquationSolverOperation::SolveEquations); - if (allocatedAuxMemory) { - this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations); - } + uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(std::vector& x, std::vector const& b) const { // Get a Jacobi decomposition of the matrix A (if not already available). - if(jacobiDecomposition == nullptr) { - std::pair, std::vector> nativeJacobiDecomposition = A.getJacobiDecomposition(); + if(!jacobiDecomposition) { + std::pair, std::vector> nativeJacobiDecomposition = A->getJacobiDecomposition(); // Convert the LU matrix to gmm++'s format. jacobiDecomposition = std::make_unique, std::vector>>(*storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(nativeJacobiDecomposition.first)), std::move(nativeJacobiDecomposition.second)); } + gmm::csr_matrix const& jacobiLU = jacobiDecomposition->first; + std::vector const& jacobiD = jacobiDecomposition->second; std::vector* currentX = &x; - std::vector* nextX = auxiliaryJacobiMemory.get(); + if(!this->auxiliaryRowVector) { + this->auxiliaryRowVector = std::make_unique>(getMatrixRowCount()); + } + std::vector* nextX = this->auxiliaryRowVector.get(); + // Set up additional environment variables. uint_fast64_t iterationCount = 0; bool converged = false; while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) { // Compute D^-1 * (b - LU * x) and store result in nextX. - gmm::mult_add(jacobiDecomposition->first, gmm::scaled(*currentX, -storm::utility::one()), b, *nextX); - storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition->second, *nextX, *nextX); + gmm::mult_add(jacobiLU, gmm::scaled(*currentX, -storm::utility::one()), b, *nextX); + storm::utility::vector::multiplyVectorsPointwise(jacobiD, *nextX, *nextX); // Now check if the process already converged within our precision. converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion()); @@ -246,21 +264,16 @@ namespace storm { // 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 == auxiliaryJacobiMemory.get()) { + if (currentX == this->auxiliaryRowVector.get()) { std::swap(x, *currentX); } - // If we allocated auxiliary memory, we need to dispose of it now. - if (allocatedAuxMemory) { - this->deallocateAuxMemory(LinearEquationSolverOperation::SolveEquations); - } - return iterationCount; } template - GmmxxLinearEquationSolverSettings& GmmxxLinearEquationSolver::getSettings() { - return settings; + void GmmxxLinearEquationSolver::setSettings(GmmxxLinearEquationSolverSettings const& newSettings) { + settings = newSettings; } template @@ -269,54 +282,12 @@ namespace storm { } template - bool GmmxxLinearEquationSolver::allocateAuxMemory(LinearEquationSolverOperation operation) const { - bool result = false; - if (operation == LinearEquationSolverOperation::SolveEquations) { - if (this->getSettings().getSolutionMethod() == GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi) { - if (!auxiliaryJacobiMemory) { - auxiliaryJacobiMemory = std::make_unique>(this->getMatrixRowCount()); - result = true; - } - } - } - result |= LinearEquationSolver::allocateAuxMemory(operation); - return result; - } - - template - bool GmmxxLinearEquationSolver::deallocateAuxMemory(LinearEquationSolverOperation operation) const { - bool result = false; - if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliaryJacobiMemory) { - result = true; - auxiliaryJacobiMemory.reset(); - } - } - result |= LinearEquationSolver::deallocateAuxMemory(operation); - return result; - } - - template - bool GmmxxLinearEquationSolver::reallocateAuxMemory(LinearEquationSolverOperation operation) const { - bool result = false; - if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliaryJacobiMemory) { - result = auxiliaryJacobiMemory->size() != this->getMatrixColumnCount(); - auxiliaryJacobiMemory->resize(this->getMatrixRowCount()); - } - } - result |= LinearEquationSolver::reallocateAuxMemory(operation); - return result; - } - - template - bool GmmxxLinearEquationSolver::hasAuxMemory(LinearEquationSolverOperation operation) const { - bool result = false; - if (operation == LinearEquationSolverOperation::SolveEquations) { - result |= static_cast(auxiliaryJacobiMemory); - } - result |= LinearEquationSolver::hasAuxMemory(operation); - return result; + void GmmxxLinearEquationSolver::resetAuxiliaryData() const { + gmmxxA.reset(); + iluPreconditioner.reset(); + diagonalPreconditioner.reset(); + jacobiDecomposition.reset(); + LinearEquationSolver::resetAuxiliaryData(); } template diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index 825c35a3d..09a50a92e 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -95,29 +95,26 @@ namespace storm { virtual bool solveEquations(std::vector& x, std::vector const& b) const override; virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; - GmmxxLinearEquationSolverSettings& getSettings(); + void setSettings(GmmxxLinearEquationSolverSettings const& newSettings); GmmxxLinearEquationSolverSettings const& getSettings() const; - virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const override; - virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const override; - virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const override; - virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const override; + /* + * Clears auxiliary data that has possibly been stored during previous calls of the solver. + */ + virtual void resetAuxiliaryData() 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 set to zero. * @param b The right-hand side of the equation system. * @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) const; + uint_fast64_t solveLinearEquationSystemWithJacobi(std::vector& x, std::vector const& b) const; virtual uint64_t getMatrixRowCount() const override; virtual uint64_t getMatrixColumnCount() const override; @@ -130,14 +127,13 @@ namespace storm { // the pointer refers to localA. storm::storage::SparseMatrix const* A; - // The (gmm++) matrix associated with this equation solver. - std::unique_ptr> gmmxxMatrix; - // The settings used by the solver. GmmxxLinearEquationSolverSettings settings; - // Auxiliary storage for the Jacobi method. - mutable std::unique_ptr> auxiliaryJacobiMemory; + // Auxiliary data obtained during solving + mutable std::unique_ptr> gmmxxA; + mutable std::unique_ptr>> iluPreconditioner; + mutable std::unique_ptr>> diagonalPreconditioner; mutable std::unique_ptr, std::vector>> jacobiDecomposition; }; diff --git a/src/solver/LinearEquationSolver.cpp b/src/solver/LinearEquationSolver.cpp index c4d46cd46..029e596fc 100644 --- a/src/solver/LinearEquationSolver.cpp +++ b/src/solver/LinearEquationSolver.cpp @@ -14,21 +14,21 @@ namespace storm { namespace solver { template - LinearEquationSolver::LinearEquationSolver() : auxiliaryRepeatedMultiplyMemory(nullptr) { + LinearEquationSolver::LinearEquationSolver() { // Intentionally left empty. } template void LinearEquationSolver::repeatedMultiply(std::vector& x, std::vector const* b, uint_fast64_t n) const { - bool allocatedAuxMemory = !this->hasAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly); - if (allocatedAuxMemory) { - this->allocateAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly); + + if(!auxiliaryRowVector) { + auxiliaryRowVector = std::make_unique>(getMatrixRowCount()); } // Set up some temporary variables so that we can just swap pointers instead of copying the result after // each iteration. std::vector* currentX = &x; - std::vector* nextX = auxiliaryRepeatedMultiplyMemory.get(); + std::vector* nextX = auxiliaryRowVector.get(); // Now perform matrix-vector multiplication as long as we meet the bound. for (uint_fast64_t i = 0; i < n; ++i) { @@ -38,56 +38,16 @@ namespace storm { // 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 == auxiliaryRepeatedMultiplyMemory.get()) { + if (currentX == auxiliaryRowVector.get()) { std::swap(x, *currentX); } - - // If we allocated auxiliary memory, we need to dispose of it now. - if (allocatedAuxMemory) { - this->deallocateAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly); - } } template - bool LinearEquationSolver::allocateAuxMemory(LinearEquationSolverOperation operation) const { - if (!auxiliaryRepeatedMultiplyMemory) { - auxiliaryRepeatedMultiplyMemory = std::make_unique>(this->getMatrixColumnCount()); - return true; - } - return false; + void LinearEquationSolver::resetAuxiliaryData() const { + auxiliaryRowVector.reset(); } - - template - bool LinearEquationSolver::deallocateAuxMemory(LinearEquationSolverOperation operation) const { - if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { - if (auxiliaryRepeatedMultiplyMemory) { - auxiliaryRepeatedMultiplyMemory.reset(); - return true; - } - } - return false; - } - - template - bool LinearEquationSolver::reallocateAuxMemory(LinearEquationSolverOperation operation) const { - bool result = false; - if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { - if (auxiliaryRepeatedMultiplyMemory) { - result = auxiliaryRepeatedMultiplyMemory->size() != this->getMatrixColumnCount(); - auxiliaryRepeatedMultiplyMemory->resize(this->getMatrixColumnCount()); - } - } - return result; - } - - template - bool LinearEquationSolver::hasAuxMemory(LinearEquationSolverOperation operation) const { - if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { - return static_cast(auxiliaryRepeatedMultiplyMemory); - } - return false; - } - + template std::unique_ptr> LinearEquationSolverFactory::create(storm::storage::SparseMatrix&& matrix) const { return create(matrix); diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h index b52778051..7988d2c6a 100644 --- a/src/solver/LinearEquationSolver.h +++ b/src/solver/LinearEquationSolver.h @@ -68,40 +68,15 @@ namespace storm { */ void repeatedMultiply(std::vector& x, std::vector const* b, uint_fast64_t n) const; - // Methods related to allocating/freeing auxiliary storage. - - /*! - * Allocates auxiliary memory that can be used to perform the provided operation. Repeated calls to the - * corresponding function can then be run without allocating/deallocating this memory repeatedly. - * Note: Since the allocated memory is fit to the currently selected options of the solver, they must not - * be changed any more after allocating the auxiliary memory until it is deallocated again. - * - * @return True iff auxiliary memory was allocated. - */ - virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const; - - /*! - * Destroys previously allocated auxiliary memory for the provided operation. - * - * @return True iff auxiliary memory was deallocated. - */ - virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const; - - /*! - * If the matrix dimensions changed and auxiliary memory was allocated, this function needs to be called to - * update the auxiliary memory. - * - * @return True iff the auxiliary memory was reallocated. - */ - virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const; - - /*! - * Checks whether the solver has allocated auxiliary memory for the provided operation. - * - * @return True iff auxiliary memory was previously allocated (and not yet deallocated). + /* + * Clears auxiliary data that has possibly been stored during previous calls of the solver. */ - virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const; + virtual void resetAuxiliaryData() const; + protected: + // auxiliary storage. If set, this vector has getMatrixRowCount() entries. + mutable std::unique_ptr> auxiliaryRowVector; + private: /*! * Retrieves the row count of the matrix associated with this solver. @@ -113,8 +88,6 @@ namespace storm { */ virtual uint64_t getMatrixColumnCount() const = 0; - // Auxiliary memory for repeated matrix-vector multiplication. - mutable std::unique_ptr> auxiliaryRepeatedMultiplyMemory; }; template diff --git a/src/solver/MinMaxLinearEquationSolver.h b/src/solver/MinMaxLinearEquationSolver.h index c17aa269c..c563d8d36 100644 --- a/src/solver/MinMaxLinearEquationSolver.h +++ b/src/solver/MinMaxLinearEquationSolver.h @@ -129,7 +129,7 @@ namespace storm { virtual bool getRelative() const = 0; /* - * Resets the auxiliary data that has been stored during previous calls of this solver + * Clears auxiliary data that has possibly been stored during previous calls of the solver. */ virtual void resetAuxiliaryData() const; diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp index 42e861f21..31fe74227 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -84,12 +84,12 @@ namespace storm { } template - NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingMemory(nullptr) { + NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings) { this->setMatrix(A); } template - NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingMemory(nullptr) { + NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings) { this->setMatrix(std::move(A)); } @@ -97,20 +97,21 @@ namespace storm { void NativeLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { localA.reset(); this->A = &A; + resetAuxiliaryData(); } template void NativeLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { localA = std::make_unique>(std::move(A)); this->A = localA.get(); + resetAuxiliaryData(); } template bool NativeLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { - bool allocatedAuxStorage = !this->hasAuxMemory(LinearEquationSolverOperation::SolveEquations); - if (allocatedAuxStorage) { - this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations); + if(!this->auxiliaryRowVector) { + this->auxiliaryRowVector = std::make_unique>(getMatrixRowCount()); } if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::GaussSeidel) { @@ -125,29 +126,29 @@ namespace storm { A->performSuccessiveOverRelaxationStep(omega, x, b); // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision(*auxiliarySolvingMemory, x, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)); + converged = storm::utility::vector::equalModuloPrecision(*this->auxiliaryRowVector, x, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)); // If we did not yet converge, we need to backup the contents of x. if (!converged) { - *auxiliarySolvingMemory = x; + *this->auxiliaryRowVector = x; } // Increase iteration count so we can abort if convergence is too slow. ++iterationCount; } - // If we allocated auxiliary memory, we need to dispose of it now. - if (allocatedAuxStorage) { - this->deallocateAuxMemory(LinearEquationSolverOperation::SolveEquations); - } return converged; } else { // Get a Jacobi decomposition of the matrix A. - std::pair, std::vector> jacobiDecomposition = A->getJacobiDecomposition(); + if(!jacobiDecomposition) { + jacobiDecomposition = std::make_unique, std::vector>>(A->getJacobiDecomposition()); + } + storm::storage::SparseMatrix const& jacobiLU = jacobiDecomposition->first; + std::vector const& jacobiD = jacobiDecomposition->second; std::vector* currentX = &x; - std::vector* nextX = auxiliarySolvingMemory.get(); + std::vector* nextX = this->auxiliaryRowVector.get(); // Set up additional environment variables. uint_fast64_t iterationCount = 0; @@ -155,9 +156,9 @@ namespace storm { while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) { // Compute D^-1 * (b - LU * x) and store result in nextX. - jacobiDecomposition.first.multiplyWithVector(*currentX, *nextX); + jacobiLU.multiplyWithVector(*currentX, *nextX); storm::utility::vector::subtractVectors(b, *nextX, *nextX); - storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, *nextX, *nextX); + storm::utility::vector::multiplyVectorsPointwise(jacobiD, *nextX, *nextX); // Now check if the process already converged within our precision. converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()); @@ -171,18 +172,12 @@ namespace storm { // 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 == auxiliarySolvingMemory.get()) { + if (currentX == this->auxiliaryRowVector.get()) { std::swap(x, *currentX); } - - // If we allocated auxiliary memory, we need to dispose of it now. - if (allocatedAuxStorage) { - this->deallocateAuxMemory(LinearEquationSolverOperation::SolveEquations); - } + return iterationCount < this->getSettings().getMaximalNumberOfIterations(); } - - } template @@ -194,17 +189,21 @@ namespace storm { } } else { // If the two vectors are aliases, we need to create a temporary. - std::vector tmp(result.size()); - A->multiplyWithVector(x, tmp); + if(!this->auxiliaryRowVector) { + this->auxiliaryRowVector = std::make_unique>(getMatrixRowCount()); + } + A->multiplyWithVector(x, *this->auxiliaryRowVector); if (b != nullptr) { - storm::utility::vector::addVectors(tmp, *b, result); + storm::utility::vector::addVectors(*this->auxiliaryRowVector, *b, result); + } else { + result.swap(*this->auxiliaryRowVector); } } } template - NativeLinearEquationSolverSettings& NativeLinearEquationSolver::getSettings() { - return settings; + void NativeLinearEquationSolver::setSettings(NativeLinearEquationSolverSettings const& newSettings) { + settings = newSettings; } template @@ -213,52 +212,9 @@ namespace storm { } template - bool NativeLinearEquationSolver::allocateAuxMemory(LinearEquationSolverOperation operation) const { - bool result = false; - if (operation == LinearEquationSolverOperation::SolveEquations) { - if (!auxiliarySolvingMemory) { - auxiliarySolvingMemory = std::make_unique>(this->getMatrixRowCount()); - result = true; - } - } - result |= LinearEquationSolver::allocateAuxMemory(operation); - return result; - } - - template - bool NativeLinearEquationSolver::deallocateAuxMemory(LinearEquationSolverOperation operation) const { - bool result = false; - if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliarySolvingMemory) { - result = true; - auxiliarySolvingMemory.reset(); - } - } - result |= LinearEquationSolver::deallocateAuxMemory(operation); - return result; - } - - template - bool NativeLinearEquationSolver::reallocateAuxMemory(LinearEquationSolverOperation operation) const { - bool result = false; - if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliarySolvingMemory) { - result = auxiliarySolvingMemory->size() != this->getMatrixColumnCount(); - auxiliarySolvingMemory->resize(this->getMatrixRowCount()); - } - } - result |= LinearEquationSolver::reallocateAuxMemory(operation); - return result; - } - - template - bool NativeLinearEquationSolver::hasAuxMemory(LinearEquationSolverOperation operation) const { - bool result = false; - if (operation == LinearEquationSolverOperation::SolveEquations) { - result |= static_cast(auxiliarySolvingMemory); - } - result |= LinearEquationSolver::hasAuxMemory(operation); - return result; + void NativeLinearEquationSolver::resetAuxiliaryData() const { + jacobiDecomposition.reset(); + LinearEquationSolver::resetAuxiliaryData(); } template @@ -301,4 +257,4 @@ namespace storm { template class NativeLinearEquationSolver; template class NativeLinearEquationSolverFactory; } -} \ No newline at end of file +} diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index 4bf9275dd..df1b91a13 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -52,15 +52,14 @@ namespace storm { virtual bool solveEquations(std::vector& x, std::vector const& b) const override; virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; - NativeLinearEquationSolverSettings& getSettings(); + void setSettings(NativeLinearEquationSolverSettings const& newSettings); NativeLinearEquationSolverSettings const& getSettings() const; - virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const override; - virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const override; - virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const override; - virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const override; + /* + * Clears auxiliary data that has possibly been stored during previous calls of the solver. + */ + virtual void resetAuxiliaryData() const override; - storm::storage::SparseMatrix const& getMatrix() const { return *A; } private: virtual uint64_t getMatrixRowCount() const override; virtual uint64_t getMatrixColumnCount() const override; @@ -76,8 +75,7 @@ namespace storm { // The settings used by the solver. NativeLinearEquationSolverSettings settings; - // Auxiliary memory for the equation solving methods. - mutable std::unique_ptr> auxiliarySolvingMemory; + mutable std::unique_ptr, std::vector>> jacobiDecomposition; }; template diff --git a/src/solver/StandardMinMaxLinearEquationSolver.cpp b/src/solver/StandardMinMaxLinearEquationSolver.cpp index 360146bf2..454b59b3d 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/solver/StandardMinMaxLinearEquationSolver.cpp @@ -99,10 +99,10 @@ namespace storm { std::vector scheduler(this->A.getRowGroupCount()); // Get a vector for storing the right-hand side of the inner equation system. - if(!auxiliaryData.rowGroupVector) { - auxiliaryData.rowGroupVector = std::make_unique>(this->A.getRowGroupCount()); + if(!auxiliaryRowGroupVector) { + auxiliaryRowGroupVector = std::make_unique>(this->A.getRowGroupCount()); } - std::vector& subB = *auxiliaryData.rowGroupVector; + std::vector& subB = *auxiliaryRowGroupVector; // Resolve the nondeterminism according to the current scheduler. storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(scheduler, true); @@ -111,7 +111,6 @@ namespace storm { // Create a solver that we will use throughout the procedure. We will modify the matrix in each iteration. auto solver = linearEquationSolverFactory->create(std::move(submatrix)); - solver->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations); Status status = Status::InProgress; uint64_t iterations = 0; @@ -139,6 +138,8 @@ namespace storm { choiceValue += b[choice]; // If the value is strictly better than the solution of the inner system, we need to improve the scheduler. + // TODO: If the underlying solver is not precise, this might run forever (i.e. when a state has two choices where the (exact) values are equal). + // only changing the scheduler if the values are not equal (modulo precision) would make this unsound. if (valueImproved(dir, x[group], choiceValue)) { schedulerImproved = true; scheduler[group] = choice - this->A.getRowGroupIndices()[group]; @@ -205,19 +206,19 @@ namespace storm { template bool StandardMinMaxLinearEquationSolver::solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { - if(!auxiliaryData.linEqSolver) { - auxiliaryData.linEqSolver = linearEquationSolverFactory->create(A); + if(!linEqSolverA) { + linEqSolverA = linearEquationSolverFactory->create(A); } - if (!auxiliaryData.rowVector.get()) { - auxiliaryData.rowVector = std::make_unique>(A.getRowCount()); + if (!auxiliaryRowVector.get()) { + auxiliaryRowVector = std::make_unique>(A.getRowCount()); } - std::vector& multiplyResult = *auxiliaryData.rowVector; + std::vector& multiplyResult = *auxiliaryRowVector; - if (!auxiliaryData.rowGroupVector.get()) { - auxiliaryData.rowGroupVector = std::make_unique>(A.getRowGroupCount()); + if (!auxiliaryRowGroupVector.get()) { + auxiliaryRowGroupVector = std::make_unique>(A.getRowGroupCount()); } - std::vector* newX = auxiliaryData.rowGroupVector.get(); + std::vector* newX = auxiliaryRowGroupVector.get(); std::vector* currentX = &x; @@ -227,7 +228,7 @@ namespace storm { Status status = Status::InProgress; while (status == Status::InProgress) { // Compute x' = A*x + b. - auxiliaryData.linEqSolver->multiply(*currentX, &b, multiplyResult); + linEqSolverA->multiply(*currentX, &b, multiplyResult); // Reduce the vector x' by applying min/max for all non-deterministic choices. storm::utility::vector::reduceVectorMinOrMax(dir, multiplyResult, *newX, this->A.getRowGroupIndices()); @@ -248,14 +249,14 @@ namespace storm { // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result // is currently stored in currentX, but x is the output vector. - if (currentX == auxiliaryData.rowGroupVector.get()) { + if (currentX == auxiliaryRowGroupVector.get()) { std::swap(x, *currentX); } // If requested, we store the scheduler for retrieval. if (this->isTrackSchedulerSet()) { if(iterations==0){ //may happen due to custom termination condition. Then we need to compute x'= A*x+b - auxiliaryData.linEqSolver->multiply(x, &b, multiplyResult); + linEqSolverA->multiply(x, &b, multiplyResult); } std::vector choices(this->A.getRowGroupCount()); // Reduce the multiplyResult and keep track of the choices made @@ -272,17 +273,17 @@ namespace storm { template void StandardMinMaxLinearEquationSolver::repeatedMultiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const { - if(!auxiliaryData.linEqSolver) { - auxiliaryData.linEqSolver = linearEquationSolverFactory->create(A); + if(!linEqSolverA) { + linEqSolverA = linearEquationSolverFactory->create(A); } - if (!auxiliaryData.rowVector.get()) { - auxiliaryData.rowVector = std::make_unique>(A.getRowCount()); + if (!auxiliaryRowVector.get()) { + auxiliaryRowVector = std::make_unique>(A.getRowCount()); } - std::vector& multiplyResult = *auxiliaryData.rowVector; + std::vector& multiplyResult = *auxiliaryRowVector; for (uint64_t i = 0; i < n; ++i) { - auxiliaryData.linEqSolver->multiply(x, b, multiplyResult); + linEqSolverA->multiply(x, b, multiplyResult); // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost // element of the min/max operator stack. @@ -325,9 +326,10 @@ namespace storm { template void StandardMinMaxLinearEquationSolver::resetAuxiliaryData() const { - auxiliaryData.linEqSolver.reset(); - auxiliaryData.rowVector.reset(); - auxiliaryData.rowGroupVector.reset(); + linEqSolverA.reset(); + auxiliaryRowVector.reset(); + auxiliaryRowGroupVector.reset(); + MinMaxLinearEquationSolver::resetAuxiliaryData(); } template diff --git a/src/solver/StandardMinMaxLinearEquationSolver.h b/src/solver/StandardMinMaxLinearEquationSolver.h index d460887be..dac6a45e0 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/solver/StandardMinMaxLinearEquationSolver.h @@ -60,13 +60,9 @@ namespace storm { Converged, TerminatedEarly, MaximalIterationsExceeded, InProgress }; - struct AuxiliaryData { - std::unique_ptr> linEqSolver; - std::unique_ptr> rowVector; // A.rowCount() entries - std::unique_ptr> rowGroupVector; // A.rowGroupCount() entries - }; - - mutable AuxiliaryData auxiliaryData; + mutable std::unique_ptr> linEqSolverA; + mutable std::unique_ptr> auxiliaryRowVector; // A.rowCount() entries + mutable std::unique_ptr> auxiliaryRowGroupVector; // A.rowGroupCount() entries Status updateStatusIfNotConverged(Status status, std::vector const& x, uint64_t iterations) const; void reportStatus(Status status, uint64_t iterations) const; diff --git a/src/utility/policyguessing.cpp b/src/utility/policyguessing.cpp index 072662c33..15335b36a 100644 --- a/src/utility/policyguessing.cpp +++ b/src/utility/policyguessing.cpp @@ -171,14 +171,16 @@ namespace storm { storm::utility::vector::selectVectorValues(subB, probGreater0States, b); std::unique_ptr> linEqSysSolver(static_cast*>(storm::solver::GmmxxLinearEquationSolverFactory().create(eqSysA).release())); - auto& eqSettings = linEqSysSolver->getSettings(); + auto eqSettings = linEqSysSolver->getSettings(); eqSettings.setRelativeTerminationCriterion(relative); eqSettings.setMaximalNumberOfIterations(500); + linEqSysSolver->setSettings(eqSettings); std::size_t iterations = 0; std::vector copyX(subX.size()); ValueType precisionChangeFactor = storm::utility::one(); do { eqSettings.setPrecision(eqSettings.getPrecision() * precisionChangeFactor); + linEqSysSolver->setSettings(eqSettings); if(!linEqSysSolver->solveEquations(subX, subB)){ // break; //Solver did not converge.. so we have to go on with the current solution. } diff --git a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp index 4cf07d6ec..ccb9f8cbb 100644 --- a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp +++ b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp @@ -54,11 +54,13 @@ TEST(GmmxxLinearEquationSolver, gmres) { std::vector b = {16, -4, -7}; storm::solver::GmmxxLinearEquationSolver solver(A); - solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Gmres); - solver.getSettings().setPrecision(1e-6); - solver.getSettings().setMaximalNumberOfIterations(10000); - solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); - solver.getSettings().setNumberOfIterationsUntilRestart(50); + auto settings = solver.getSettings(); + settings.setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Gmres); + settings.setPrecision(1e-6); + settings.setMaximalNumberOfIterations(10000); + settings.setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + settings.setNumberOfIterationsUntilRestart(50); + solver.setSettings(settings); ASSERT_NO_THROW(solver.solveEquations(x, b)); ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); @@ -86,16 +88,20 @@ TEST(GmmxxLinearEquationSolver, qmr) { std::vector b = {16, -4, -7}; storm::solver::GmmxxLinearEquationSolver solver(A); - solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Qmr); - solver.getSettings().setPrecision(1e-6); - solver.getSettings().setMaximalNumberOfIterations(10000); - solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + auto settings = solver.getSettings(); + settings.setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Qmr); + settings.setPrecision(1e-6); + settings.setMaximalNumberOfIterations(10000); + settings.setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + solver.setSettings(settings); storm::solver::GmmxxLinearEquationSolver solver2(A); - solver2.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Qmr); - solver2.getSettings().setPrecision(1e-6); - solver2.getSettings().setMaximalNumberOfIterations(10000); - solver2.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + auto settings2 = solver2.getSettings(); + settings2.setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Qmr); + settings2.setPrecision(1e-6); + settings2.setMaximalNumberOfIterations(10000); + settings2.setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + solver2.setSettings(settings2); ASSERT_NO_THROW(solver2.solveEquations(x, b)); ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); @@ -123,10 +129,12 @@ TEST(GmmxxLinearEquationSolver, bicgstab) { std::vector b = {16, -4, -7}; storm::solver::GmmxxLinearEquationSolver solver(A); - solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Bicgstab); - solver.getSettings().setPrecision(1e-6); - solver.getSettings().setMaximalNumberOfIterations(10000); - solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + auto settings = solver.getSettings(); + settings.setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Bicgstab); + settings.setPrecision(1e-6); + settings.setMaximalNumberOfIterations(10000); + settings.setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + solver.setSettings(settings); ASSERT_NO_THROW(solver.solveEquations(x, b)); ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); @@ -154,10 +162,12 @@ TEST(GmmxxLinearEquationSolver, jacobi) { std::vector b = {11, -16, 1}; storm::solver::GmmxxLinearEquationSolver solver(A); - solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi); - solver.getSettings().setPrecision(1e-6); - solver.getSettings().setMaximalNumberOfIterations(10000); - solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + auto settings = solver.getSettings(); + settings.setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi); + settings.setPrecision(1e-6); + settings.setMaximalNumberOfIterations(10000); + settings.setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + solver.setSettings(settings); ASSERT_NO_THROW(solver.solveEquations(x, b)); ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); @@ -184,11 +194,13 @@ TEST(GmmxxLinearEquationSolver, gmresilu) { std::vector b = {16, -4, -7}; storm::solver::GmmxxLinearEquationSolver solver(A); - solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Gmres); - solver.getSettings().setPrecision(1e-6); - solver.getSettings().setMaximalNumberOfIterations(10000); - solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::Ilu); - solver.getSettings().setNumberOfIterationsUntilRestart(50); + auto settings = solver.getSettings(); + settings.setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Gmres); + settings.setPrecision(1e-6); + settings.setMaximalNumberOfIterations(10000); + settings.setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::Ilu); + settings.setNumberOfIterationsUntilRestart(50); + solver.setSettings(settings); ASSERT_NO_THROW(solver.solveEquations(x, b)); ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); @@ -215,11 +227,13 @@ TEST(GmmxxLinearEquationSolver, gmresdiag) { std::vector b = {16, -4, -7}; storm::solver::GmmxxLinearEquationSolver solver(A); - solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Gmres); - solver.getSettings().setPrecision(1e-6); - solver.getSettings().setMaximalNumberOfIterations(10000); - solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::Diagonal); - solver.getSettings().setNumberOfIterationsUntilRestart(50); + auto settings = solver.getSettings(); + settings.setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Gmres); + settings.setPrecision(1e-6); + settings.setMaximalNumberOfIterations(10000); + settings.setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::Diagonal); + settings.setNumberOfIterationsUntilRestart(50); + solver.setSettings(settings); ASSERT_NO_THROW(solver.solveEquations(x, b)); ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); @@ -248,4 +262,4 @@ TEST(GmmxxLinearEquationSolver, MatrixVectorMultiplication) { storm::solver::GmmxxLinearEquationSolver solver(A); ASSERT_NO_THROW(solver.repeatedMultiply(x, nullptr, 4)); ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); -} \ No newline at end of file +}