dehnert
12 years ago
15 changed files with 1270 additions and 1247 deletions
-
2CMakeLists.txt
-
309src/modelchecker/prctl/GmmxxDtmcPrctlModelChecker.h
-
186src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h
-
69src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
-
1319src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
-
17src/models/AtomicPropositionsLabeling.h
-
26src/solver/AbstractLinearEquationSolver.h
-
175src/solver/AbstractNondeterministicLinearEquationSolver.h
-
228src/solver/GmmxxLinearEquationSolver.h
-
122src/solver/GmmxxNondeterministicLinearEquationSolver.h
-
20src/storage/BitVector.h
-
9src/storm.cpp
-
10src/utility/Settings.cpp
-
10src/utility/Settings.h
-
17src/utility/vector.h
@ -1,309 +0,0 @@ |
|||||
/* |
|
||||
* GmmxxDtmcPrctlModelChecker.h |
|
||||
* |
|
||||
* Created on: 06.12.2012 |
|
||||
* Author: Christian Dehnert |
|
||||
*/ |
|
||||
|
|
||||
#ifndef STORM_MODELCHECKER_PRCTL_GMMXXDTMCPRCTLMODELCHECKER_H_ |
|
||||
#define STORM_MODELCHECKER_PRCTL_GMMXXDTMCPRCTLMODELCHECKER_H_ |
|
||||
|
|
||||
#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" |
|
||||
#include "src/adapters/GmmxxAdapter.h" |
|
||||
#include "src/storage/SparseMatrix.h" |
|
||||
#include "src/utility/ConstTemplates.h" |
|
||||
#include "src/utility/Settings.h" |
|
||||
|
|
||||
#include "gmm/gmm_matrix.h" |
|
||||
#include "gmm/gmm_iter_solvers.h" |
|
||||
|
|
||||
#include <cmath> |
|
||||
|
|
||||
namespace storm { |
|
||||
namespace modelchecker { |
|
||||
namespace prctl { |
|
||||
|
|
||||
/* |
|
||||
* An implementation of the SparseDtmcPrctlModelChecker interface that uses gmm++ as the solving backend. |
|
||||
*/ |
|
||||
template <class Type> |
|
||||
class GmmxxDtmcPrctlModelChecker : public SparseDtmcPrctlModelChecker<Type> { |
|
||||
|
|
||||
public: |
|
||||
/*! |
|
||||
* Constructs a GmmxxDtmcPrctlModelChecker with the given model. |
|
||||
* |
|
||||
* @param model The DTMC to be checked. |
|
||||
*/ |
|
||||
explicit GmmxxDtmcPrctlModelChecker(storm::models::Dtmc<Type> const& dtmc) : SparseDtmcPrctlModelChecker<Type>(dtmc) { |
|
||||
// Intentionally left empty. |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Copy constructs a GmmxxDtmcPrctlModelChecker from the given model checker. In particular, this means that the newly |
|
||||
* constructed model checker will have the model of the given model checker as its associated model. |
|
||||
*/ |
|
||||
explicit GmmxxDtmcPrctlModelChecker(storm::modelchecker::prctl::GmmxxDtmcPrctlModelChecker<Type> const& modelchecker) : SparseDtmcPrctlModelChecker<Type>(modelchecker) { |
|
||||
// Intentionally left empty. |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Virtual destructor. Needs to be virtual, because this class has virtual methods. |
|
||||
*/ |
|
||||
virtual ~GmmxxDtmcPrctlModelChecker() { |
|
||||
// Intentionally left empty. |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Returns the name of this module. |
|
||||
* @returns The name of this module. |
|
||||
*/ |
|
||||
static std::string getModuleName() { |
|
||||
return "gmm++det"; |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Returns a trigger such that if the option "matrixlib" is set to "gmm++", this model checker |
|
||||
* is to be used. |
|
||||
* @returns An option trigger for this module. |
|
||||
*/ |
|
||||
static std::pair<std::string, std::string> getOptionTrigger() { |
|
||||
return std::pair<std::string, std::string>("matrixlib", "gmm++"); |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Registers all options associated with the gmm++ matrix library. |
|
||||
*/ |
|
||||
static void putOptions(boost::program_options::options_description* desc) { |
|
||||
desc->add_options()("lemethod", boost::program_options::value<std::string>()->default_value("bicgstab")->notifier(&validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr, jacobi}."); |
|
||||
desc->add_options()("maxiter", boost::program_options::value<unsigned>()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving."); |
|
||||
desc->add_options()("precision", boost::program_options::value<double>()->default_value(1e-6), "Sets the precision for iterative equation solving."); |
|
||||
desc->add_options()("precond", boost::program_options::value<std::string>()->default_value("ilu")->notifier(&validatePreconditioner), "Sets the preconditioning technique for linear equation solving. Must be in {ilu, diagonal, ildlt, none}."); |
|
||||
desc->add_options()("relative", boost::program_options::value<bool>()->default_value(true), "Sets whether the relative or absolute error is considered for detecting convergence."); |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Validates whether the given lemethod matches one of the available ones. |
|
||||
* Throws an exception of type InvalidSettings in case the selected method is illegal. |
|
||||
*/ |
|
||||
static void validateLeMethod(const std::string& lemethod) { |
|
||||
if ((lemethod != "bicgstab") && (lemethod != "qmr") && (lemethod != "jacobi")) { |
|
||||
throw exceptions::InvalidSettingsException() << "Argument " << lemethod << " for option 'lemethod' is invalid."; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Validates whether the given preconditioner matches one of the available ones. |
|
||||
* Throws an exception of type InvalidSettings in case the selected preconditioner is illegal. |
|
||||
*/ |
|
||||
static void validatePreconditioner(const std::string& preconditioner) { |
|
||||
if ((preconditioner != "ilu") && (preconditioner != "diagonal") && (preconditioner != "ildlt") && (preconditioner != "none")) { |
|
||||
throw exceptions::InvalidSettingsException() << "Argument " << preconditioner << " for option 'precond' is invalid."; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
private: |
|
||||
|
|
||||
virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b, uint_fast64_t n = 1) const { |
|
||||
// Transform the transition probability A to the gmm++ format to use its arithmetic. |
|
||||
gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
|
||||
|
|
||||
// Set up some temporary variables so that we can just swap pointers instead of copying the result after each |
|
||||
// iteration. |
|
||||
std::vector<Type>* swap = nullptr; |
|
||||
std::vector<Type>* currentVector = &x; |
|
||||
std::vector<Type>* tmpVector = new std::vector<Type>(this->getModel().getNumberOfStates()); |
|
||||
|
|
||||
// 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; |
|
||||
} |
|
||||
|
|
||||
virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const { |
|
||||
// Get the settings object to customize linear solving. |
|
||||
storm::settings::Settings* s = storm::settings::instance(); |
|
||||
|
|
||||
// Prepare an iteration object that determines the accuracy, maximum number of iterations |
|
||||
// and the like. |
|
||||
gmm::iteration iter(s->get<double>("precision"), 0, s->get<unsigned>("maxiter")); |
|
||||
|
|
||||
// Print some information about the used preconditioner. |
|
||||
const std::string& precond = s->getString("precond"); |
|
||||
LOG4CPLUS_INFO(logger, "Starting iterative solver."); |
|
||||
if (s->getString("lemethod") == "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 (s->getString("lemethod") == "bicgstab") { |
|
||||
LOG4CPLUS_INFO(logger, "Using BiCGStab method."); |
|
||||
// Transform the transition probability matrix to the gmm++ format to use its arithmetic. |
|
||||
gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
|
||||
if (precond == "ilu") { |
|
||||
gmm::bicgstab(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), iter); |
|
||||
} else if (precond == "diagonal") { |
|
||||
gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), iter); |
|
||||
} else if (precond == "ildlt") { |
|
||||
gmm::bicgstab(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*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 (s->getString("lemethod") == "qmr") { |
|
||||
LOG4CPLUS_INFO(logger, "Using QMR method."); |
|
||||
// Transform the transition probability matrix to the gmm++ format to use its arithmetic. |
|
||||
gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
|
||||
if (precond == "ilu") { |
|
||||
gmm::qmr(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), iter); |
|
||||
} else if (precond == "diagonal") { |
|
||||
gmm::qmr(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), iter); |
|
||||
} else if (precond == "ildlt") { |
|
||||
gmm::qmr(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*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 (s->getString("lemethod") == "jacobi") { |
|
||||
LOG4CPLUS_INFO(logger, "Using Jacobi method."); |
|
||||
uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b); |
|
||||
|
|
||||
// Check if the solver converged and issue a warning otherwise. |
|
||||
if (iterations < s->get<unsigned>("maxiter")) { |
|
||||
LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); |
|
||||
} else { |
|
||||
LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* 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 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. |
|
||||
*/ |
|
||||
uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const { |
|
||||
// Get the settings object to customize linear solving. |
|
||||
storm::settings::Settings* s = storm::settings::instance(); |
|
||||
|
|
||||
double precision = s->get<double>("precision"); |
|
||||
uint_fast64_t maxIterations = s->get<unsigned>("maxiter"); |
|
||||
bool relative = s->get<bool>("relative"); |
|
||||
|
|
||||
// Get a Jacobi decomposition of the matrix A. |
|
||||
storm::storage::SparseMatrix<Type>::SparseJacobiDecomposition_t jacobiDecomposition = A.getJacobiDecomposition(); |
|
||||
|
|
||||
// Convert the (inverted) diagonal matrix to gmm++'s format. |
|
||||
gmm::csr_matrix<Type>* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(std::move(jacobiDecomposition.second)); |
|
||||
// Convert the LU matrix to gmm++'s format. |
|
||||
gmm::csr_matrix<Type>* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(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<Type>* xNext = new std::vector<Type>(x.size()); |
|
||||
const std::vector<Type>* xCopy = xNext; |
|
||||
std::vector<Type>* xCurrent = &x; |
|
||||
|
|
||||
// Target vector for precision calculation. |
|
||||
std::vector<Type>* residuum = new std::vector<Type>(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, *xNext); |
|
||||
// b - R * x_k (xNext contains R * x_k) -> xNext |
|
||||
gmm::add(b, gmm::scaled(*xNext, -1.0), *xNext); |
|
||||
// D^-1 * (b - R * x_k) -> xNext |
|
||||
gmm::mult(*gmmxxDiagonalInverted, *xNext, *xNext); |
|
||||
|
|
||||
// Swap xNext with xCurrent so that the next iteration can use xCurrent again without having to copy the |
|
||||
// vector. |
|
||||
std::vector<Type> *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 residuum; |
|
||||
delete gmmxxDiagonalInverted; |
|
||||
delete gmmxxLU; |
|
||||
|
|
||||
return iterationCount; |
|
||||
} |
|
||||
}; |
|
||||
|
|
||||
} // namespace prctl |
|
||||
} // namespace modelchecker |
|
||||
} // namespace storm |
|
||||
|
|
||||
#endif /* STORM_MODELCHECKER_PRCTL_GMMXXDTMCPRCTLMODELCHECKER_H_ */ |
|
||||
@ -1,186 +0,0 @@ |
|||||
/* |
|
||||
* GmmxxDtmcPrctlModelChecker.h |
|
||||
* |
|
||||
* Created on: 06.12.2012 |
|
||||
* Author: Christian Dehnert |
|
||||
*/ |
|
||||
|
|
||||
#ifndef STORM_MODELCHECKER_PRCTL_GMMXXMDPPRCTLMODELCHECKER_H_ |
|
||||
#define STORM_MODELCHECKER_PRCTL_GMMXXMDPPRCTLMODELCHECKER_H_ |
|
||||
|
|
||||
#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h" |
|
||||
#include "src/adapters/GmmxxAdapter.h" |
|
||||
|
|
||||
#include "gmm/gmm_matrix.h" |
|
||||
#include "gmm/gmm_iter_solvers.h" |
|
||||
|
|
||||
#include <cmath> |
|
||||
|
|
||||
namespace storm { |
|
||||
namespace modelchecker { |
|
||||
namespace prctl { |
|
||||
|
|
||||
/* |
|
||||
* An implementation of the SparseMdpPrctlModelChecker interface that uses gmm++ as the solving backend. |
|
||||
*/ |
|
||||
template <class Type> |
|
||||
class GmmxxMdpPrctlModelChecker : public SparseMdpPrctlModelChecker<Type> { |
|
||||
|
|
||||
public: |
|
||||
/*! |
|
||||
* Constructs a GmmxxMdpPrctlModelChecker with the given model. |
|
||||
* |
|
||||
* @param model The MDP to be checked. |
|
||||
*/ |
|
||||
explicit GmmxxMdpPrctlModelChecker(storm::models::Mdp<Type> const& model) : SparseMdpPrctlModelChecker<Type>(model) { |
|
||||
// Intentionally left empty. |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Copy constructs a GmmxxDtmcPrctlModelChecker from the given model checker. In particular, this means that the newly |
|
||||
* constructed model checker will have the model of the given model checker as its associated model. |
|
||||
*/ |
|
||||
explicit GmmxxMdpPrctlModelChecker(storm::modelchecker::prctl::GmmxxMdpPrctlModelChecker<Type> const& modelchecker) : SparseMdpPrctlModelChecker<Type>(modelchecker) { |
|
||||
// Intentionally left empty. |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Virtual destructor. Needs to be virtual, because this class has virtual methods. |
|
||||
*/ |
|
||||
virtual ~GmmxxMdpPrctlModelChecker() { |
|
||||
// Intentionally left empty. |
|
||||
} |
|
||||
|
|
||||
private: |
|
||||
|
|
||||
/*! |
|
||||
* Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes x[i+1] = A*x[i] + b |
|
||||
* until x[n], where x[0] = x. |
|
||||
* |
|
||||
* @param A The matrix that is to be multiplied against the vector. |
|
||||
* @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter, |
|
||||
* i.e. after the method returns, this vector will contain the computed values. |
|
||||
* @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix. |
|
||||
* @param b If not null, this vector is being added to the result after each matrix-vector multiplication. |
|
||||
* @param n Specifies the number of iterations the matrix-vector multiplication is performed. |
|
||||
* @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector. |
|
||||
*/ |
|
||||
virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const { |
|
||||
// Transform the transition probability matrix to the gmm++ format to use its arithmetic. |
|
||||
gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
|
||||
|
|
||||
// Create vector for result of multiplication, which is reduced to the result vector after |
|
||||
// each multiplication. |
|
||||
std::vector<Type> multiplyResult(A.getRowCount()); |
|
||||
|
|
||||
// Now perform matrix-vector multiplication as long as we meet the bound of the formula. |
|
||||
for (uint_fast64_t i = 0; i < n; ++i) { |
|
||||
gmm::mult(*gmmxxMatrix, x, multiplyResult); |
|
||||
|
|
||||
if (b != nullptr) { |
|
||||
gmm::add(*b, multiplyResult); |
|
||||
} |
|
||||
|
|
||||
if (this->minimumOperatorStack.top()) { |
|
||||
storm::utility::vector::reduceVectorMin(multiplyResult, x, nondeterministicChoiceIndices); |
|
||||
} else { |
|
||||
storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices); |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
// Delete intermediate results and return result. |
|
||||
delete gmmxxMatrix; |
|
||||
} |
|
||||
|
|
||||
/*! |
|
||||
* Solves the given equation system under the given parameters using the power method. |
|
||||
* |
|
||||
* @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise. |
|
||||
* @param A The matrix A specifying the coefficients of the equations. |
|
||||
* @param x The vector x for which to solve the equations. The initial value of the elements of |
|
||||
* this vector are used as the initial guess and might thus influence performance and convergence. |
|
||||
* @param b The vector b specifying the values on the right-hand-sides of the equations. |
|
||||
* @return The solution of the system of linear equations in form of the elements of the vector |
|
||||
* x. |
|
||||
*/ |
|
||||
void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<uint_fast64_t>* takenChoices = nullptr) const override { |
|
||||
// Get the settings object to customize solving. |
|
||||
storm::settings::Settings* s = storm::settings::instance(); |
|
||||
|
|
||||
// Get relevant user-defined settings for solving the equations. |
|
||||
double precision = s->get<double>("precision"); |
|
||||
unsigned maxIterations = s->get<unsigned>("maxiter"); |
|
||||
bool relative = s->get<bool>("relative"); |
|
||||
|
|
||||
// Transform the transition probability matrix to the gmm++ format to use its arithmetic. |
|
||||
gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
|
||||
|
|
||||
// Set up the environment for the power method. |
|
||||
std::vector<Type> multiplyResult(A.getRowCount()); |
|
||||
std::vector<Type>* currentX = &x; |
|
||||
std::vector<Type>* newX = new std::vector<Type>(x.size()); |
|
||||
std::vector<Type>* swap = nullptr; |
|
||||
uint_fast64_t iterations = 0; |
|
||||
bool converged = false; |
|
||||
|
|
||||
// Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number |
|
||||
// of iterations. |
|
||||
while (!converged && iterations < maxIterations) { |
|
||||
std::cout << "iter: " << iterations << std::endl; |
|
||||
// Compute x' = A*x + b. |
|
||||
gmm::mult(*gmmxxMatrix, *currentX, multiplyResult); |
|
||||
gmm::add(b, multiplyResult); |
|
||||
|
|
||||
// Reduce the vector x' by applying min/max for all non-deterministic choices. |
|
||||
if (minimize) { |
|
||||
storm::utility::vector::reduceVectorMin(multiplyResult, *newX, nondeterministicChoiceIndices); |
|
||||
} else { |
|
||||
storm::utility::vector::reduceVectorMax(multiplyResult, *newX, nondeterministicChoiceIndices); |
|
||||
} |
|
||||
|
|
||||
// Determine whether the method converged. |
|
||||
converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative); |
|
||||
|
|
||||
if (!converged) { |
|
||||
for (uint_fast64_t i = 0; i < newX->size(); ++i) { |
|
||||
std::cout << (*currentX)[i] << " vs. " << (*newX)[i] << std::endl; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
|
|
||||
// Update environment variables. |
|
||||
swap = currentX; |
|
||||
currentX = newX; |
|
||||
newX = swap; |
|
||||
++iterations; |
|
||||
} |
|
||||
|
|
||||
// If we were requested to record the taken choices, we have to construct the vector now. |
|
||||
if (takenChoices != nullptr) { |
|
||||
this->computeTakenChoices(minimize, multiplyResult, *takenChoices, nondeterministicChoiceIndices); |
|
||||
} |
|
||||
|
|
||||
// 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 (iterations % 2 == 1) { |
|
||||
std::swap(x, *currentX); |
|
||||
delete currentX; |
|
||||
} else { |
|
||||
delete newX; |
|
||||
} |
|
||||
delete gmmxxMatrix; |
|
||||
|
|
||||
// Check if the solver converged and issue a warning otherwise. |
|
||||
if (converged) { |
|
||||
LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); |
|
||||
} else { |
|
||||
LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations."); |
|
||||
} |
|
||||
} |
|
||||
}; |
|
||||
|
|
||||
} //namespace prctl |
|
||||
} //namespace modelchecker |
|
||||
} //namespace storm |
|
||||
|
|
||||
#endif /* STORM_MODELCHECKER_PRCTL_GMMXXMDPPRCTLMODELCHECKER_H_ */ |
|
||||
1319
src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
File diff suppressed because it is too large
View File
File diff suppressed because it is too large
View File
@ -0,0 +1,26 @@ |
|||||
|
#ifndef STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_ |
||||
|
#define STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_ |
||||
|
|
||||
|
#include "src/storage/SparseMatrix.h" |
||||
|
|
||||
|
#include <vector> |
||||
|
|
||||
|
namespace storm { |
||||
|
namespace solver { |
||||
|
|
||||
|
template<class Type> |
||||
|
class AbstractLinearEquationSolver { |
||||
|
public: |
||||
|
|
||||
|
virtual AbstractLinearEquationSolver<Type>* clone() const = 0; |
||||
|
|
||||
|
virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const = 0; |
||||
|
|
||||
|
virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const = 0; |
||||
|
|
||||
|
}; |
||||
|
|
||||
|
} // namespace solver |
||||
|
} // namespace storm |
||||
|
|
||||
|
#endif /* STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_ */ |
||||
@ -0,0 +1,175 @@ |
|||||
|
#ifndef STORM_SOLVER_ABSTRACTNONDETERMINISTICLINEAREQUATIONSOLVER_H_ |
||||
|
#define STORM_SOLVER_ABSTRACTNONDETERMINISTICLINEAREQUATIONSOLVER_H_ |
||||
|
|
||||
|
#include "src/storage/SparseMatrix.h" |
||||
|
|
||||
|
#include <vector> |
||||
|
|
||||
|
namespace storm { |
||||
|
namespace solver { |
||||
|
|
||||
|
template<class Type> |
||||
|
class AbstractNondeterministicLinearEquationSolver { |
||||
|
public: |
||||
|
|
||||
|
virtual AbstractNondeterministicLinearEquationSolver<Type>* clone() const { |
||||
|
return new AbstractNondeterministicLinearEquationSolver<Type>(); |
||||
|
} |
||||
|
|
||||
|
/*! |
||||
|
* Solves the equation system A*x = b given by the parameters. |
||||
|
* |
||||
|
* @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise. |
||||
|
* @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 b The right-hand side of the equation system. |
||||
|
* @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix. |
||||
|
* @param takenChoices If not null, this vector will be filled with the nondeterministic choices taken by the states |
||||
|
* to achieve the solution of the equation system. This assumes that the given vector has at least as many elements |
||||
|
* as there are states in the MDP. |
||||
|
* @returns The solution vector x of the system of linear equations as the content of the parameter x. |
||||
|
*/ |
||||
|
virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const { |
||||
|
LOG4CPLUS_INFO(logger, "Starting iterative solver."); |
||||
|
|
||||
|
// Get the settings object to customize solving. |
||||
|
storm::settings::Settings* s = storm::settings::instance(); |
||||
|
|
||||
|
// Get relevant user-defined settings for solving the equations. |
||||
|
double precision = s->get<double>("precision"); |
||||
|
unsigned maxIterations = s->get<unsigned>("maxiter"); |
||||
|
bool relative = s->get<bool>("relative"); |
||||
|
|
||||
|
// Set up the environment for the power method. |
||||
|
std::vector<Type> multiplyResult(A.getRowCount()); |
||||
|
std::vector<Type>* currentX = &x; |
||||
|
std::vector<Type>* newX = new std::vector<Type>(x.size()); |
||||
|
std::vector<Type>* swap = nullptr; |
||||
|
uint_fast64_t iterations = 0; |
||||
|
bool converged = false; |
||||
|
|
||||
|
// Proceed with the iterations as long as the method did not converge or reach the |
||||
|
// user-specified maximum number of iterations. |
||||
|
while (!converged && iterations < maxIterations) { |
||||
|
// Compute x' = A*x + b. |
||||
|
A.multiplyWithVector(*currentX, multiplyResult); |
||||
|
storm::utility::vector::addVectorsInPlace(multiplyResult, b); |
||||
|
|
||||
|
// 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. |
||||
|
if (minimize) { |
||||
|
storm::utility::vector::reduceVectorMin(multiplyResult, *newX, nondeterministicChoiceIndices); |
||||
|
} else { |
||||
|
storm::utility::vector::reduceVectorMax(multiplyResult, *newX, nondeterministicChoiceIndices); |
||||
|
} |
||||
|
|
||||
|
// Determine whether the method converged. |
||||
|
converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative); |
||||
|
|
||||
|
// Update environment variables. |
||||
|
swap = currentX; |
||||
|
currentX = newX; |
||||
|
newX = swap; |
||||
|
++iterations; |
||||
|
} |
||||
|
|
||||
|
// 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 (iterations % 2 == 1) { |
||||
|
std::swap(x, *currentX); |
||||
|
delete currentX; |
||||
|
} else { |
||||
|
delete newX; |
||||
|
} |
||||
|
|
||||
|
// Check if the solver converged and issue a warning otherwise. |
||||
|
if (converged) { |
||||
|
LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); |
||||
|
} else { |
||||
|
LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/*! |
||||
|
* Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes x[i+1] = A*x[i] + b |
||||
|
* until x[n], where x[0] = x. |
||||
|
* |
||||
|
* @param minimize If set, all choices are resolved such that the solution value becomes minimal and maximal otherwise. |
||||
|
* @param A The matrix that is to be multiplied against the vector. |
||||
|
* @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter, |
||||
|
* i.e. after the method returns, this vector will contain the computed values. |
||||
|
* @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix. |
||||
|
* @param b If not null, this vector is being added to the result after each matrix-vector multiplication. |
||||
|
* @param n Specifies the number of iterations the matrix-vector multiplication is performed. |
||||
|
* @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector. |
||||
|
*/ |
||||
|
virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const { |
||||
|
// Create vector for result of multiplication, which is reduced to the result vector after |
||||
|
// each multiplication. |
||||
|
std::vector<Type> multiplyResult(A.getRowCount()); |
||||
|
|
||||
|
// Now perform matrix-vector multiplication as long as we meet the bound of the formula. |
||||
|
for (uint_fast64_t i = 0; i < n; ++i) { |
||||
|
A.multiplyWithVector(x, multiplyResult); |
||||
|
|
||||
|
// Add b if it is non-null. |
||||
|
if (b != nullptr) { |
||||
|
storm::utility::vector::addVectorsInPlace(multiplyResult, *b); |
||||
|
} |
||||
|
|
||||
|
// 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. |
||||
|
if (minimize) { |
||||
|
storm::utility::vector::reduceVectorMin(multiplyResult, x, nondeterministicChoiceIndices); |
||||
|
} else { |
||||
|
storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/*! |
||||
|
* Returns the name of this solver. |
||||
|
* |
||||
|
* @returns The name of this solver. |
||||
|
*/ |
||||
|
static std::string getName() { |
||||
|
return "native"; |
||||
|
} |
||||
|
|
||||
|
/*! |
||||
|
* Registers all options associated with the gmm++ matrix library. |
||||
|
*/ |
||||
|
static void putOptions(boost::program_options::options_description* desc) { |
||||
|
desc->add_options()("lemethod", boost::program_options::value<std::string>()->default_value("bicgstab")->notifier(&validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr, jacobi}."); |
||||
|
desc->add_options()("maxiter", boost::program_options::value<unsigned>()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving."); |
||||
|
desc->add_options()("precision", boost::program_options::value<double>()->default_value(1e-6), "Sets the precision for iterative equation solving."); |
||||
|
desc->add_options()("precond", boost::program_options::value<std::string>()->default_value("ilu")->notifier(&validatePreconditioner), "Sets the preconditioning technique for linear equation solving. Must be in {ilu, diagonal, ildlt, none}."); |
||||
|
desc->add_options()("relative", boost::program_options::value<bool>()->default_value(true), "Sets whether the relative or absolute error is considered for detecting convergence."); |
||||
|
} |
||||
|
|
||||
|
/*! |
||||
|
* Validates whether the given lemethod matches one of the available ones. |
||||
|
* Throws an exception of type InvalidSettings in case the selected method is illegal. |
||||
|
*/ |
||||
|
static void validateLeMethod(const std::string& lemethod) { |
||||
|
if ((lemethod != "bicgstab") && (lemethod != "qmr") && (lemethod != "jacobi")) { |
||||
|
throw exceptions::InvalidSettingsException() << "Argument " << lemethod << " for option 'lemethod' is invalid."; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/*! |
||||
|
* Validates whether the given preconditioner matches one of the available ones. |
||||
|
* Throws an exception of type InvalidSettings in case the selected preconditioner is illegal. |
||||
|
*/ |
||||
|
static void validatePreconditioner(const std::string& preconditioner) { |
||||
|
if ((preconditioner != "ilu") && (preconditioner != "diagonal") && (preconditioner != "ildlt") && (preconditioner != "none")) { |
||||
|
throw exceptions::InvalidSettingsException() << "Argument " << preconditioner << " for option 'precond' is invalid."; |
||||
|
} |
||||
|
} |
||||
|
}; |
||||
|
|
||||
|
} // namespace solver |
||||
|
} // namespace storm |
||||
|
|
||||
|
#endif /* STORM_SOLVER_ABSTRACTNONDETERMINISTICLINEAREQUATIONSOLVER_H_ */ |
||||
@ -0,0 +1,228 @@ |
|||||
|
#ifndef STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_ |
||||
|
#define STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_ |
||||
|
|
||||
|
#include "AbstractLinearEquationSolver.h" |
||||
|
#include "src/adapters/GmmxxAdapter.h" |
||||
|
#include "src/utility/ConstTemplates.h" |
||||
|
#include "src/utility/Settings.h" |
||||
|
|
||||
|
#include "gmm/gmm_matrix.h" |
||||
|
#include "gmm/gmm_iter_solvers.h" |
||||
|
|
||||
|
#include <cmath> |
||||
|
|
||||
|
namespace storm { |
||||
|
namespace solver { |
||||
|
|
||||
|
template<class Type> |
||||
|
class GmmxxLinearEquationSolver : public AbstractLinearEquationSolver<Type> { |
||||
|
public: |
||||
|
|
||||
|
virtual AbstractLinearEquationSolver<Type>* clone() const { |
||||
|
return new GmmxxLinearEquationSolver<Type>(); |
||||
|
} |
||||
|
|
||||
|
virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const { |
||||
|
// Get the settings object to customize linear solving. |
||||
|
storm::settings::Settings* s = storm::settings::instance(); |
||||
|
|
||||
|
// Prepare an iteration object that determines the accuracy, maximum number of iterations |
||||
|
// and the like. |
||||
|
gmm::iteration iter(s->get<double>("precision"), 0, s->get<unsigned>("maxiter")); |
||||
|
|
||||
|
// Print some information about the used preconditioner. |
||||
|
const std::string& precond = s->getString("precond"); |
||||
|
LOG4CPLUS_INFO(logger, "Starting iterative solver."); |
||||
|
if (s->getString("lemethod") == "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 (s->getString("lemethod") == "bicgstab") { |
||||
|
LOG4CPLUS_INFO(logger, "Using BiCGStab method."); |
||||
|
// Transform the transition probability matrix to the gmm++ format to use its arithmetic. |
||||
|
gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
||||
|
if (precond == "ilu") { |
||||
|
gmm::bicgstab(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), iter); |
||||
|
} else if (precond == "diagonal") { |
||||
|
gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), iter); |
||||
|
} else if (precond == "ildlt") { |
||||
|
gmm::bicgstab(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*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 (s->getString("lemethod") == "qmr") { |
||||
|
LOG4CPLUS_INFO(logger, "Using QMR method."); |
||||
|
// Transform the transition probability matrix to the gmm++ format to use its arithmetic. |
||||
|
gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
||||
|
if (precond == "ilu") { |
||||
|
gmm::qmr(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), iter); |
||||
|
} else if (precond == "diagonal") { |
||||
|
gmm::qmr(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), iter); |
||||
|
} else if (precond == "ildlt") { |
||||
|
gmm::qmr(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*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 (s->getString("lemethod") == "jacobi") { |
||||
|
LOG4CPLUS_INFO(logger, "Using Jacobi method."); |
||||
|
uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b); |
||||
|
|
||||
|
// Check if the solver converged and issue a warning otherwise. |
||||
|
if (iterations < s->get<unsigned>("maxiter")) { |
||||
|
LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); |
||||
|
} else { |
||||
|
LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b, uint_fast64_t n = 1) const { |
||||
|
// Transform the transition probability A to the gmm++ format to use its arithmetic. |
||||
|
gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
||||
|
|
||||
|
// Set up some temporary variables so that we can just swap pointers instead of copying the result after each |
||||
|
// iteration. |
||||
|
std::vector<Type>* swap = nullptr; |
||||
|
std::vector<Type>* currentVector = &x; |
||||
|
std::vector<Type>* tmpVector = new std::vector<Type>(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; |
||||
|
} |
||||
|
|
||||
|
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 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. |
||||
|
*/ |
||||
|
uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const { |
||||
|
// Get the settings object to customize linear solving. |
||||
|
storm::settings::Settings* s = storm::settings::instance(); |
||||
|
|
||||
|
double precision = s->get<double>("precision"); |
||||
|
uint_fast64_t maxIterations = s->get<unsigned>("maxiter"); |
||||
|
bool relative = s->get<bool>("relative"); |
||||
|
|
||||
|
// Get a Jacobi decomposition of the matrix A. |
||||
|
typename storm::storage::SparseMatrix<Type>::SparseJacobiDecomposition_t jacobiDecomposition = A.getJacobiDecomposition(); |
||||
|
|
||||
|
// Convert the (inverted) diagonal matrix to gmm++'s format. |
||||
|
gmm::csr_matrix<Type>* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(std::move(jacobiDecomposition.second)); |
||||
|
// Convert the LU matrix to gmm++'s format. |
||||
|
gmm::csr_matrix<Type>* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(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<Type>* xNext = new std::vector<Type>(x.size()); |
||||
|
const std::vector<Type>* xCopy = xNext; |
||||
|
std::vector<Type>* xCurrent = &x; |
||||
|
|
||||
|
// Target vector for precision calculation. |
||||
|
std::vector<Type>* residuum = new std::vector<Type>(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, *xNext); |
||||
|
// b - R * x_k (xNext contains R * x_k) -> xNext |
||||
|
gmm::add(b, gmm::scaled(*xNext, -1.0), *xNext); |
||||
|
// D^-1 * (b - R * x_k) -> xNext |
||||
|
gmm::mult(*gmmxxDiagonalInverted, *xNext, *xNext); |
||||
|
|
||||
|
// Swap xNext with xCurrent so that the next iteration can use xCurrent again without having to copy the |
||||
|
// vector. |
||||
|
std::vector<Type> *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 residuum; |
||||
|
delete gmmxxDiagonalInverted; |
||||
|
delete gmmxxLU; |
||||
|
|
||||
|
return iterationCount; |
||||
|
} |
||||
|
|
||||
|
}; |
||||
|
|
||||
|
} // namespace solver |
||||
|
} // namespace storm |
||||
|
|
||||
|
#endif /* STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_ */ |
||||
@ -0,0 +1,122 @@ |
|||||
|
#ifndef STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_ |
||||
|
#define STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_ |
||||
|
|
||||
|
#include "AbstractLinearEquationSolver.h" |
||||
|
#include "src/adapters/GmmxxAdapter.h" |
||||
|
#include "src/storage/SparseMatrix.h" |
||||
|
|
||||
|
#include "gmm/gmm_matrix.h" |
||||
|
|
||||
|
#include <vector> |
||||
|
|
||||
|
namespace storm { |
||||
|
namespace solver { |
||||
|
|
||||
|
template<class Type> |
||||
|
class GmmxxNondeterministicLinearEquationSolver : public AbstractNondeterministicLinearEquationSolver<Type> { |
||||
|
public: |
||||
|
|
||||
|
virtual AbstractLinearEquationSolver<Type>* clone() const { |
||||
|
return new GmmxxNondeterministicLinearEquationSolver<Type>(); |
||||
|
} |
||||
|
|
||||
|
virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const { |
||||
|
// Transform the transition probability matrix to the gmm++ format to use its arithmetic. |
||||
|
gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
||||
|
|
||||
|
// Create vector for result of multiplication, which is reduced to the result vector after |
||||
|
// each multiplication. |
||||
|
std::vector<Type> multiplyResult(A.getRowCount()); |
||||
|
|
||||
|
// Now perform matrix-vector multiplication as long as we meet the bound of the formula. |
||||
|
for (uint_fast64_t i = 0; i < n; ++i) { |
||||
|
gmm::mult(*gmmxxMatrix, x, multiplyResult); |
||||
|
|
||||
|
if (b != nullptr) { |
||||
|
gmm::add(*b, multiplyResult); |
||||
|
} |
||||
|
|
||||
|
if (this->minimumOperatorStack.top()) { |
||||
|
storm::utility::vector::reduceVectorMin(multiplyResult, x, nondeterministicChoiceIndices); |
||||
|
} else { |
||||
|
storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Delete intermediate results and return result. |
||||
|
delete gmmxxMatrix; |
||||
|
} |
||||
|
|
||||
|
void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<uint_fast64_t>* takenChoices = nullptr) const override { |
||||
|
// Get the settings object to customize solving. |
||||
|
storm::settings::Settings* s = storm::settings::instance(); |
||||
|
|
||||
|
// Get relevant user-defined settings for solving the equations. |
||||
|
double precision = s->get<double>("precision"); |
||||
|
unsigned maxIterations = s->get<unsigned>("maxiter"); |
||||
|
bool relative = s->get<bool>("relative"); |
||||
|
|
||||
|
// Transform the transition probability matrix to the gmm++ format to use its arithmetic. |
||||
|
gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A); |
||||
|
|
||||
|
// Set up the environment for the power method. |
||||
|
std::vector<Type> multiplyResult(A.getRowCount()); |
||||
|
std::vector<Type>* currentX = &x; |
||||
|
std::vector<Type>* newX = new std::vector<Type>(x.size()); |
||||
|
std::vector<Type>* swap = nullptr; |
||||
|
uint_fast64_t iterations = 0; |
||||
|
bool converged = false; |
||||
|
|
||||
|
// Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number |
||||
|
// of iterations. |
||||
|
while (!converged && iterations < maxIterations) { |
||||
|
// Compute x' = A*x + b. |
||||
|
gmm::mult(*gmmxxMatrix, *currentX, multiplyResult); |
||||
|
gmm::add(b, multiplyResult); |
||||
|
|
||||
|
// Reduce the vector x' by applying min/max for all non-deterministic choices. |
||||
|
if (minimize) { |
||||
|
storm::utility::vector::reduceVectorMin(multiplyResult, *newX, nondeterministicChoiceIndices); |
||||
|
} else { |
||||
|
storm::utility::vector::reduceVectorMax(multiplyResult, *newX, nondeterministicChoiceIndices); |
||||
|
} |
||||
|
|
||||
|
// Determine whether the method converged. |
||||
|
converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative); |
||||
|
|
||||
|
// Update environment variables. |
||||
|
swap = currentX; |
||||
|
currentX = newX; |
||||
|
newX = swap; |
||||
|
++iterations; |
||||
|
} |
||||
|
|
||||
|
// If we were requested to record the taken choices, we have to construct the vector now. |
||||
|
if (takenChoices != nullptr) { |
||||
|
this->computeTakenChoices(minimize, multiplyResult, *takenChoices, nondeterministicChoiceIndices); |
||||
|
} |
||||
|
|
||||
|
// 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 (iterations % 2 == 1) { |
||||
|
std::swap(x, *currentX); |
||||
|
delete currentX; |
||||
|
} else { |
||||
|
delete newX; |
||||
|
} |
||||
|
delete gmmxxMatrix; |
||||
|
|
||||
|
// Check if the solver converged and issue a warning otherwise. |
||||
|
if (converged) { |
||||
|
LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); |
||||
|
} else { |
||||
|
LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations."); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
}; |
||||
|
|
||||
|
} // namespace solver |
||||
|
} // namespace storm |
||||
|
|
||||
|
#endif /* STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_ */ |
||||
Write
Preview
Loading…
Cancel
Save
Reference in new issue