From 2e8f74a70090ac7754de88845bb762a6158abf80 Mon Sep 17 00:00:00 2001 From: dehnert Date: Mon, 1 Feb 2016 14:51:34 +0100 Subject: [PATCH] First steps toward exact linear equation solver Former-commit-id: 669af13b84574d35c6df40e16ca8a007ed303dc4 --- src/adapters/CarlAdapter.h | 9 +- .../EliminationLinearEquationSolver.cpp | 62 +++++++++++++ src/solver/EliminationLinearEquationSolver.h | 34 +++++++ src/solver/NativeLinearEquationSolver.h | 2 +- src/storage/FlexibleSparseMatrix.cpp | 88 +++++++++++++++++++ src/storage/FlexibleSparseMatrix.h | 65 ++++++++++++++ .../expressions/ToRationalFunctionVisitor.cpp | 4 +- src/utility/vector.h | 2 +- 8 files changed, 260 insertions(+), 6 deletions(-) create mode 100644 src/solver/EliminationLinearEquationSolver.cpp create mode 100644 src/solver/EliminationLinearEquationSolver.h create mode 100644 src/storage/FlexibleSparseMatrix.cpp create mode 100644 src/storage/FlexibleSparseMatrix.h diff --git a/src/adapters/CarlAdapter.h b/src/adapters/CarlAdapter.h index 2ed172cde..093ea6ae7 100644 --- a/src/adapters/CarlAdapter.h +++ b/src/adapters/CarlAdapter.h @@ -4,6 +4,8 @@ // Include config to know whether CARL is available or not. #include "storm-config.h" +#include + #ifdef STORM_HAVE_CARL #include @@ -42,9 +44,12 @@ namespace carl { } namespace storm { - typedef cln::cl_RA RationalNumber; +// typedef boost::multiprecision::gmp_rational RationalNumber; + typedef mpq_class RationalNumber; + + typedef cln::cl_RA CarlRationalNumber; typedef carl::Variable Variable; - typedef carl::MultivariatePolynomial RawPolynomial; + typedef carl::MultivariatePolynomial RawPolynomial; typedef carl::FactorizedPolynomial Polynomial; typedef carl::Relation CompareRelation; diff --git a/src/solver/EliminationLinearEquationSolver.cpp b/src/solver/EliminationLinearEquationSolver.cpp new file mode 100644 index 000000000..67a281b70 --- /dev/null +++ b/src/solver/EliminationLinearEquationSolver.cpp @@ -0,0 +1,62 @@ +#include "src/solver/EliminationLinearEquationSolver.h" + +#include "src/utility/vector.h" + +namespace storm { + namespace solver { + template + EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix const& A) : A(A) { + // Intentionally left empty. + } + + template + void EliminationLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + // TODO: implement state-elimination here. + } + + template + void EliminationLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { + // Set up some temporary variables so that we can just swap pointers instead of copying the result after + // each iteration. + std::vector* currentX = &x; + + bool multiplyResultProvided = true; + std::vector* nextX = multiplyResult; + if (nextX == nullptr) { + nextX = new std::vector(x.size()); + multiplyResultProvided = false; + } + std::vector const* copyX = nextX; + + // Now perform matrix-vector multiplication as long as we meet the bound. + for (uint_fast64_t i = 0; i < n; ++i) { + A.multiplyWithVector(*currentX, *nextX); + std::swap(nextX, currentX); + + // If requested, add an offset to the current result vector. + if (b != nullptr) { + storm::utility::vector::addVectors(*currentX, *b, *currentX); + } + } + + // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, + // because the output is supposed to be stored in the input vector x. + if (currentX == copyX) { + std::swap(x, *currentX); + } + + // If the vector for the temporary multiplication result was not provided, we need to delete it. + if (!multiplyResultProvided) { + delete copyX; + } + } + + template class EliminationLinearEquationSolver; + + // TODO: make this work with the proper implementation of solveEquationSystem. + template class EliminationLinearEquationSolver; + template class EliminationLinearEquationSolver; + + } +} + diff --git a/src/solver/EliminationLinearEquationSolver.h b/src/solver/EliminationLinearEquationSolver.h new file mode 100644 index 000000000..25dffa651 --- /dev/null +++ b/src/solver/EliminationLinearEquationSolver.h @@ -0,0 +1,34 @@ +#ifndef STORM_SOLVER_ELIMINATIONLINEAREQUATIONSOLVER_H_ +#define STORM_SOLVER_ELIMINATIONLINEAREQUATIONSOLVER_H_ + +#include "src/solver/LinearEquationSolver.h" + +namespace storm { + namespace solver { + /*! + * A class that uses gaussian elimination to implement the LinearEquationSolver interface. In particular + */ + template + class EliminationLinearEquationSolver : public LinearEquationSolver { + public: + + /*! + * Constructs a linear equation solver. + * + * @param A The matrix defining the coefficients of the linear equation system. + */ + EliminationLinearEquationSolver(storm::storage::SparseMatrix const& A); + + virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; + + virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; + + private: + // A reference to the original matrix used for this equation solver. + storm::storage::SparseMatrix const& A; + + }; + } +} + +#endif /* STORM_SOLVER_ELIMINATIONLINEAREQUATIONSOLVER_H_ */ \ No newline at end of file diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index 6fa972f05..42dbbd3d9 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -10,6 +10,7 @@ namespace storm { enum class NativeLinearEquationSolverSolutionMethod { Jacobi, GaussSeidel, SOR }; + /*! * A class that uses StoRM's native matrix operations to implement the LinearEquationSolver interface. */ @@ -17,7 +18,6 @@ namespace storm { class NativeLinearEquationSolver : public LinearEquationSolver { public: - /*! * Constructs a linear equation solver with parameters being set according to the settings object. * diff --git a/src/storage/FlexibleSparseMatrix.cpp b/src/storage/FlexibleSparseMatrix.cpp new file mode 100644 index 000000000..b77f59736 --- /dev/null +++ b/src/storage/FlexibleSparseMatrix.cpp @@ -0,0 +1,88 @@ +#include "src/storage/FlexibleSparseMatrix.h" + +#include "src/storage/SparseMatrix.h" +#include "src/storage/BitVector.h" + +namespace storm { + namespace storage { + template + FlexibleSparseMatrix::FlexibleSparseMatrix(index_type rows) : data(rows) { + // Intentionally left empty. + } + + template + void FlexibleSparseMatrix::reserveInRow(index_type row, index_type numberOfElements) { + this->data[row].reserve(numberOfElements); + } + + template + typename FlexibleSparseMatrix::row_type& FlexibleSparseMatrix::getRow(index_type index) { + return this->data[index]; + } + + template + typename FlexibleSparseMatrix::row_type const& FlexibleSparseMatrix::getRow(index_type index) const { + return this->data[index]; + } + + template + typename FlexibleSparseMatrix::index_type FlexibleSparseMatrix::getNumberOfRows() const { + return this->data.size(); + } + + template + bool FlexibleSparseMatrix::hasSelfLoop(storm::storage::sparse::state_type state) { + for (auto const& entry : this->getRow(state)) { + if (entry.getColumn() < state) { + continue; + } else if (entry.getColumn() > state) { + return false; + } else if (entry.getColumn() == state) { + return true; + } + } + return false; + } + + template + void FlexibleSparseMatrix::print() const { + for (uint_fast64_t index = 0; index < this->data.size(); ++index) { + std::cout << index << " - "; + for (auto const& element : this->getRow(index)) { + std::cout << "(" << element.getColumn() << ", " << element.getValue() << ") "; + } + std::cout << std::endl; + } + } + + template + bool FlexibleSparseMatrix::empty() const { + for (auto const& row : this->data) { + if (!row.empty()) { + return false; + } + } + return true; + } + + template + void FlexibleSparseMatrix::filter(storm::storage::BitVector const& rowFilter, storm::storage::BitVector const& columnFilter) { + for (uint_fast64_t rowIndex = 0; rowIndex < this->data.size(); ++rowIndex) { + auto& row = this->data[rowIndex]; + if (!rowFilter.get(rowIndex)) { + row.clear(); + row.shrink_to_fit(); + continue; + } + row_type newRow; + for (auto const& element : row) { + if (columnFilter.get(element.getColumn())) { + newRow.push_back(element); + } + } + row = std::move(newRow); + } + } + + } +} \ No newline at end of file diff --git a/src/storage/FlexibleSparseMatrix.h b/src/storage/FlexibleSparseMatrix.h new file mode 100644 index 000000000..6a30109ac --- /dev/null +++ b/src/storage/FlexibleSparseMatrix.h @@ -0,0 +1,65 @@ +#ifndef STORM_STORAGE_FLEXIBLESPARSEMATRIX_H_ +#define STORM_STORAGE_FLEXIBLESPARSEMATRIX_H_ + +#include +#include + +#include "src/storage/sparse/StateType.h" + +namespace storm { + namespace storage { + template + class MatrixEntry; + + class BitVector; + + template + class FlexibleSparseMatrix { + public: + // TODO: make this class a bit more consistent with the big sparse matrix and improve it: + // * rename getNumberOfRows -> getRowCount + // * store number of columns to also provide getColumnCount + // * rename hasSelfLoop -> rowHasDiagonalElement + // * add output iterator and improve the way the matrix is printed + // * add conversion functionality from/to sparse matrix + // * add documentation + // * rename filter to something more appropriate (getSubmatrix?) + // * add stuff like clearRow, multiplyRowWithScalar + + typedef uint_fast64_t index_type; + typedef ValueType value_type; + typedef std::vector> row_type; + typedef typename row_type::iterator iterator; + typedef typename row_type::const_iterator const_iterator; + + FlexibleSparseMatrix() = default; + FlexibleSparseMatrix(index_type rows); + + void reserveInRow(index_type row, index_type numberOfElements); + + row_type& getRow(index_type); + row_type const& getRow(index_type) const; + + index_type getNumberOfRows() const; + + void print() const; + + bool empty() const; + + void filter(storm::storage::BitVector const& rowFilter, storm::storage::BitVector const& columnFilter); + + /*! + * Checks whether the given state has a self-loop with an arbitrary probability in the probability matrix. + * + * @param state The state for which to check whether it possesses a self-loop. + * @return True iff the given state has a self-loop with an arbitrary probability in the probability matrix. + */ + bool hasSelfLoop(storm::storage::sparse::state_type state); + + private: + std::vector data; + }; + } +} + +#endif /* STORM_STORAGE_FLEXIBLESPARSEMATRIX_H_ */ \ No newline at end of file diff --git a/src/storage/expressions/ToRationalFunctionVisitor.cpp b/src/storage/expressions/ToRationalFunctionVisitor.cpp index 97501aa1b..5e0ebcd64 100644 --- a/src/storage/expressions/ToRationalFunctionVisitor.cpp +++ b/src/storage/expressions/ToRationalFunctionVisitor.cpp @@ -88,12 +88,12 @@ namespace storm { template boost::any ToRationalFunctionVisitor::visit(IntegerLiteralExpression const& expression) { - return RationalFunctionType(carl::rationalize(static_cast(expression.getValue()))); + return RationalFunctionType(carl::rationalize(static_cast(expression.getValue()))); } template boost::any ToRationalFunctionVisitor::visit(DoubleLiteralExpression const& expression) { - return RationalFunctionType(carl::rationalize(expression.getValue())); + return RationalFunctionType(carl::rationalize(expression.getValue())); } template class ToRationalFunctionVisitor; diff --git a/src/utility/vector.h b/src/utility/vector.h index 665569645..a4e551efa 100644 --- a/src/utility/vector.h +++ b/src/utility/vector.h @@ -269,7 +269,7 @@ namespace storm { * @param target The target vector. */ template - void applyPointwise(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, std::function function) { + void applyPointwise(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, std::function const& function) { #ifdef STORM_HAVE_INTELTBB tbb::parallel_for(tbb::blocked_range(0, target.size()), [&](tbb::blocked_range const& range) {