From 31cbe14d3c2aa834f44a1073dabf26f1dce30f1e Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Mon, 9 Mar 2020 13:01:10 +0100 Subject: [PATCH] Multiplier: Added a flag to specify whether gaussSeidel style multiplications should be performed forward or backwards. --- .../solver/AcyclicLinearEquationSolver.h | 47 +++++++++ src/storm/solver/GmmxxMultiplier.cpp | 97 ++++++++++++++----- src/storm/solver/GmmxxMultiplier.h | 8 +- src/storm/solver/Multiplier.cpp | 4 +- src/storm/solver/Multiplier.h | 8 +- src/storm/solver/NativeMultiplier.cpp | 16 ++- src/storm/solver/NativeMultiplier.h | 4 +- 7 files changed, 145 insertions(+), 39 deletions(-) create mode 100644 src/storm/solver/AcyclicLinearEquationSolver.h diff --git a/src/storm/solver/AcyclicLinearEquationSolver.h b/src/storm/solver/AcyclicLinearEquationSolver.h new file mode 100644 index 000000000..e6e5fe3a3 --- /dev/null +++ b/src/storm/solver/AcyclicLinearEquationSolver.h @@ -0,0 +1,47 @@ +#pragma once + +#include "storm/solver/StandardMinMaxLinearEquationSolver.h" + +namespace storm { + + class Environment; + + namespace solver { + + /*! + * This solver can be used on equation systems that are known to be acyclic. + * It is optimized for solving many instances of the equation system with the same underlying matrix. + */ + template + class AcyclicMinMaxLinearEquationSolver : public StandardMinMaxLinearEquationSolver { + public: + AcyclicMinMaxLinearEquationSolver(); + AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A); + AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A); + + virtual ~AcyclicMinMaxLinearEquationSolver() { + } + + virtual void clearCache() const override; + + virtual MinMaxLinearEquationSolverRequirements getRequirements(Environment const& env, boost::optional const& direction = boost::none, bool const& hasInitialScheduler = false) const override ; + + protected: + + virtual bool internalSolveEquations(storm::Environment const& env, OptimizationDirection d, std::vector& x, std::vector const& b) const override; + + private: + // cached multiplier either with original matrix or ordered matrix + mutable std::unique_ptr> multiplier; + // cached matrix for the multiplier (only if different from original matrix) + mutable boost::optional> orderedMatrix; + // cached row group ordering (only if not identity) + mutable boost::optional> rowGroupOrdering; // A.rowGroupCount() entries + // can be used if the entries in 'b' need to be reordered + mutable boost::optional> auxiliaryRowVector; // A.rowCount() entries + // contains factors applied to scale the entries of the 'b' vector + mutable std::vector> bFactors; + + }; + } +} diff --git a/src/storm/solver/GmmxxMultiplier.cpp b/src/storm/solver/GmmxxMultiplier.cpp index 59d71af8e..d8c4c2468 100644 --- a/src/storm/solver/GmmxxMultiplier.cpp +++ b/src/storm/solver/GmmxxMultiplier.cpp @@ -65,13 +65,21 @@ namespace storm { } template - void GmmxxMultiplier::multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b) const { + void GmmxxMultiplier::multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b, bool backwards) const { initialize(); STORM_LOG_ASSERT(gmmMatrix.nr == gmmMatrix.nc, "Expecting square matrix."); - if (b) { - gmm::mult_add_by_row_bwd(gmmMatrix, x, *b, x, gmm::abstract_dense()); + if (backwards) { + if (b) { + gmm::mult_add_by_row_bwd(gmmMatrix, x, *b, x, gmm::abstract_dense()); + } else { + gmm::mult_by_row_bwd(gmmMatrix, x, x, gmm::abstract_dense()); + } } else { - gmm::mult_by_row_bwd(gmmMatrix, x, x, gmm::abstract_dense()); + if (b) { + gmm::mult_add_by_row(gmmMatrix, x, *b, x, gmm::abstract_dense()); + } else { + gmm::mult_by_row(gmmMatrix, x, x, gmm::abstract_dense()); + } } } @@ -90,7 +98,7 @@ namespace storm { if (parallelize(env)) { multAddReduceParallel(dir, rowGroupIndices, x, b, *target, choices); } else { - multAddReduceHelper(dir, rowGroupIndices, x, b, *target, choices); + multAddReduceHelper(dir, rowGroupIndices, x, b, *target, choices, false); } if (&x == &result) { std::swap(result, *this->cachedVector); @@ -98,9 +106,9 @@ namespace storm { } template - void GmmxxMultiplier::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices) const { + void GmmxxMultiplier::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices, bool backwards) const { initialize(); - multAddReduceHelper(dir, rowGroupIndices, x, b, x, choices); + multAddReduceHelper(dir, rowGroupIndices, x, b, x, choices, backwards); } template @@ -119,16 +127,24 @@ namespace storm { } template - void GmmxxMultiplier::multAddReduceHelper(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { + void GmmxxMultiplier::multAddReduceHelper(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices, bool backwards) const { if (dir == storm::OptimizationDirection::Minimize) { - multAddReduceHelper>(rowGroupIndices, x, b, result, choices); + if (backwards) { + multAddReduceHelper, true>(rowGroupIndices, x, b, result, choices); + } else { + multAddReduceHelper, false>(rowGroupIndices, x, b, result, choices); + } } else { - multAddReduceHelper>(rowGroupIndices, x, b, result, choices); + if (backwards) { + multAddReduceHelper, true>(rowGroupIndices, x, b, result, choices); + } else { + multAddReduceHelper, false>(rowGroupIndices, x, b, result, choices); + } } } template - template + template void GmmxxMultiplier::multAddReduceHelper(std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { Compare compare; typedef std::vector VectorType; @@ -136,29 +152,31 @@ namespace storm { typename gmm::linalg_traits::const_iterator add_it, add_ite; if (b) { - add_it = gmm::vect_end(*b) - 1; - add_ite = gmm::vect_begin(*b) - 1; + add_it = backwards ? gmm::vect_end(*b) - 1 : gmm::vect_begin(*b); + add_ite = backwards ?gmm::vect_begin(*b) - 1 : gmm::vect_end(*b); } - typename gmm::linalg_traits::iterator target_it = gmm::vect_end(result) - 1; - typename gmm::linalg_traits::const_row_iterator itr = mat_row_const_end(gmmMatrix) - 1; + typename gmm::linalg_traits::iterator target_it = backwards ? gmm::vect_end(result) - 1 : gmm::vect_begin(result); + typename gmm::linalg_traits::const_row_iterator itr = backwards ? mat_row_const_end(gmmMatrix) - 1 : mat_row_const_begin(gmmMatrix); typename std::vector::iterator choice_it; if (choices) { - choice_it = choices->end() - 1; + choice_it = backwards ? choices->end() - 1 : choices->begin(); } // Variables for correctly tracking choices (only update if new choice is strictly better). ValueType oldSelectedChoiceValue; uint64_t selectedChoice; - uint64_t currentRow = gmmMatrix.nrows() - 1; - for (auto row_group_it = rowGroupIndices.end() - 2, row_group_ite = rowGroupIndices.begin() - 1; row_group_it != row_group_ite; --row_group_it, --choice_it, --target_it) { + uint64_t currentRow = backwards ? gmmMatrix.nrows() - 1 : 0; + auto row_group_it = backwards ? rowGroupIndices.end() - 2 : rowGroupIndices.begin(); + auto row_group_ite = backwards ? rowGroupIndices.begin() - 1 : rowGroupIndices.end() - 1; + while (row_group_it != row_group_ite) { ValueType currentValue = storm::utility::zero(); // Only multiply and reduce if the row group is not empty. if (*row_group_it != *(row_group_it + 1)) { + // Process the (backwards ? last : first) row of the current row group if (b) { currentValue = *add_it; - --add_it; } currentValue += vect_sp(gmm::linalg_traits::row(itr), x); @@ -169,11 +187,21 @@ namespace storm { oldSelectedChoiceValue = currentValue; } } - - --itr; - --currentRow; - for (uint64_t row = *row_group_it + 1, rowEnd = *(row_group_it + 1); row < rowEnd; ++row, --currentRow, --itr, --add_it) { + // move row-based iterators to the next row + if (backwards) { + --itr; + --currentRow; + --add_it; + } else { + ++itr; + ++currentRow; + ++add_it; + } + + // Process the (rowGroupSize-1) remaining rows within the current row Group + uint64_t rowGroupSize = *(row_group_it + 1) - *row_group_it; + for (uint64_t i = 1; i < rowGroupSize; ++i) { ValueType newValue = b ? *add_it : storm::utility::zero(); newValue += vect_sp(gmm::linalg_traits::row(itr), x); @@ -187,6 +215,16 @@ namespace storm { selectedChoice = currentRow - *row_group_it; } } + // move row-based iterators to the next row + if (backwards) { + --itr; + --currentRow; + --add_it; + } else { + ++itr; + ++currentRow; + ++add_it; + } } // Finally write value to target vector. @@ -195,11 +233,22 @@ namespace storm { *choice_it = selectedChoice; } } + + // move rowGroup-based iterators to the next row group + if (backwards) { + --row_group_it; + --choice_it; + --target_it; + } else { + ++row_group_it; + ++choice_it; + ++target_it; + } } } template<> - template + template void GmmxxMultiplier::multAddReduceHelper(std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Operation not supported for this data type."); } diff --git a/src/storm/solver/GmmxxMultiplier.h b/src/storm/solver/GmmxxMultiplier.h index da70e2a05..f0394770b 100644 --- a/src/storm/solver/GmmxxMultiplier.h +++ b/src/storm/solver/GmmxxMultiplier.h @@ -22,9 +22,9 @@ namespace storm { virtual ~GmmxxMultiplier() = default; virtual void multiply(Environment const& env, std::vector const& x, std::vector const* b, std::vector& result) const override; - virtual void multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b) const override; + virtual void multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b, bool backwards = true) const override; virtual void multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const override; - virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices = nullptr) const override; + virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices = nullptr, bool backwards = true) const override; virtual void multiplyRow(uint64_t const& rowIndex, std::vector const& x, ValueType& value) const override; virtual void clearCache() const override; @@ -36,9 +36,9 @@ namespace storm { void multAdd(std::vector const& x, std::vector const* b, std::vector& result) const; void multAddParallel(std::vector const& x, std::vector const* b, std::vector& result) const; void multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const; - void multAddReduceHelper(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const; + void multAddReduceHelper(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr, bool backwards = true) const; - template + template void multAddReduceHelper(std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const; mutable gmm::csr_matrix gmmMatrix; diff --git a/src/storm/solver/Multiplier.cpp b/src/storm/solver/Multiplier.cpp index abcb53f3d..f0a5c4204 100644 --- a/src/storm/solver/Multiplier.cpp +++ b/src/storm/solver/Multiplier.cpp @@ -33,8 +33,8 @@ namespace storm { } template - void Multiplier::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector& x, std::vector const* b, std::vector* choices) const { - multiplyAndReduceGaussSeidel(env, dir, this->matrix.getRowGroupIndices(), x, b, choices); + void Multiplier::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector& x, std::vector const* b, std::vector* choices, bool backwards) const { + multiplyAndReduceGaussSeidel(env, dir, this->matrix.getRowGroupIndices(), x, b, choices, backwards); } template diff --git a/src/storm/solver/Multiplier.h b/src/storm/solver/Multiplier.h index 5f9036022..58aeb4e38 100644 --- a/src/storm/solver/Multiplier.h +++ b/src/storm/solver/Multiplier.h @@ -49,8 +49,9 @@ namespace storm { * to the number of columns of A. * @param b If non-null, this vector is added after the multiplication. If given, its length must be equal * to the number of rows of A. + * @param backwards if true, the iterations will be performed beginning from the last row and ending at the first row. */ - virtual void multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b) const = 0; + virtual void multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b, bool backwards = true) const = 0; /*! * Performs a matrix-vector multiplication x' = A*x + b and then minimizes/maximizes over the row groups @@ -82,9 +83,10 @@ namespace storm { * @param result The target vector into which to write the multiplication result. Its length must be equal * to the number of rows of A. Can be the same as the x vector. * @param choices If given, the choices made in the reduction process are written to this vector. + * @param backwards if true, the iterations will be performed beginning from the last rowgroup and ending at the first rowgroup. */ - void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector& x, std::vector const* b, std::vector* choices = nullptr) const; - virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices = nullptr) const = 0; + void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector& x, std::vector const* b, std::vector* choices = nullptr, bool backwards = true) const; + virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices = nullptr, bool backwards = true) const = 0; /*! * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After diff --git a/src/storm/solver/NativeMultiplier.cpp b/src/storm/solver/NativeMultiplier.cpp index b8c6c3549..8b13d07cc 100644 --- a/src/storm/solver/NativeMultiplier.cpp +++ b/src/storm/solver/NativeMultiplier.cpp @@ -53,8 +53,12 @@ namespace storm { } template - void NativeMultiplier::multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b) const { - this->matrix.multiplyWithVectorBackward(x, x, b); + void NativeMultiplier::multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b, bool backwards) const { + if (backwards) { + this->matrix.multiplyWithVectorBackward(x, x, b); + } else { + this->matrix.multiplyWithVectorForward(x, x, b); + } } template @@ -79,8 +83,12 @@ namespace storm { } template - void NativeMultiplier::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices) const { - this->matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices); + void NativeMultiplier::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices, bool backwards) const { + if (backwards) { + this->matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices); + } else { + this->matrix.multiplyAndReduceForward(dir, rowGroupIndices, x, b, x, choices); + } } template diff --git a/src/storm/solver/NativeMultiplier.h b/src/storm/solver/NativeMultiplier.h index 3d9e31a02..92424f91c 100644 --- a/src/storm/solver/NativeMultiplier.h +++ b/src/storm/solver/NativeMultiplier.h @@ -19,9 +19,9 @@ namespace storm { virtual ~NativeMultiplier() = default; virtual void multiply(Environment const& env, std::vector const& x, std::vector const* b, std::vector& result) const override; - virtual void multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b) const override; + virtual void multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b, bool backwards = true) const override; virtual void multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const override; - virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices = nullptr) const override; + virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices = nullptr, bool backwards = true) const override; virtual void multiplyRow(uint64_t const& rowIndex, std::vector const& x, ValueType& value) const override; virtual void multiplyRow2(uint64_t const& rowIndex, std::vector const& x1, ValueType& val1, std::vector const& x2, ValueType& val2) const override;