From 50245d3d86b4786b545f33c3eeab240b7a50d5e4 Mon Sep 17 00:00:00 2001 From: TimQu Date: Wed, 28 Feb 2018 12:46:21 +0100 Subject: [PATCH] gmmm multiplier --- src/storm/solver/GmmxxMultiplier.cpp | 176 +++++++++++++++++---------- src/storm/solver/GmmxxMultiplier.h | 40 +++--- 2 files changed, 137 insertions(+), 79 deletions(-) diff --git a/src/storm/solver/GmmxxMultiplier.cpp b/src/storm/solver/GmmxxMultiplier.cpp index 5158d03fb..739cbdbe5 100644 --- a/src/storm/solver/GmmxxMultiplier.cpp +++ b/src/storm/solver/GmmxxMultiplier.cpp @@ -2,6 +2,9 @@ #include "storm/adapters/RationalNumberAdapter.h" #include "storm/adapters/RationalFunctionAdapter.h" +#include "storm/storage/SparseMatrix.h" +#include "storm/settings/SettingsManager.h" +#include "storm/settings/modules/CoreSettings.h" #include "storm/utility/constants.h" #include "storm/exceptions/NotSupportedException.h" @@ -11,60 +14,113 @@ namespace storm { namespace solver { - template - GmmxxMultiplier::GmmxxMultiplier() : storm::utility::VectorHelper() { + template + GmmxxMultiplier::GmmxxMultiplier(storm::storage::SparseMatrix const& matrix) : Multiplier(matrix) { // Intentionally left empty. } - template - void GmmxxMultiplier::multAdd(gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result) const { - if (this->parallelize()) { - multAddParallel(matrix, x, b, result); - } else { - if (b) { - gmm::mult_add(matrix, x, *b, result); + template + void GmmxxMultiplier::initialize() const { + if (gmmMatrix.nrows() == 0) { + gmmMatrix = std::move(*storm::adapters::GmmxxAdapter().toGmmxxSparseMatrix(this->matrix)); + } + } + + template + void GmmxxMultiplier::clearCache() const { + gmmMatrix = gmm::csr_matrix(); + Multiplier::clearCache(); + } + + template + bool GmmxxMultiplier::parallelize(Environment const& env) const { +#ifdef STORM_HAVE_INTELTBB + return storm::settings::getModule().isUseIntelTbbSet(); +#else + return false; +#endif + } + + template + void GmmxxMultiplier::multiply(Environment const& env, std::vector const& x, std::vector const* b, std::vector& result) const { + initialize(); + std::vector* target = &result; + if (&x == &result) { + if (this->cachedVector) { + this->cachedVector->resize(x.size()); } else { - gmm::mult(matrix, x, result); + this->cachedVector = std::make_unique>(x.size()); } + target = this->cachedVector.get(); + } + if (parallelize(env)) { + multAddParallel(x, b, *target); + } else { + multAdd(x, b, *target); + } + if (&x == &result) { + std::swap(result, *this->cachedVector); } } - template - void GmmxxMultiplier::multAddGaussSeidelBackward(gmm::csr_matrix const& matrix, std::vector& x, std::vector const* b) const { - STORM_LOG_ASSERT(matrix.nr == matrix.nc, "Expecting square matrix."); + template + void GmmxxMultiplier::multiplyGaussSeidel(Environment const& env, std::vector& x, std::vector const* b) const { + initialize(); + STORM_LOG_ASSERT(gmmMatrix.nr == gmmMatrix.nc, "Expecting square matrix."); if (b) { - gmm::mult_add_by_row_bwd(matrix, x, *b, x, gmm::abstract_dense()); + gmm::mult_add_by_row_bwd(gmmMatrix, x, *b, x, gmm::abstract_dense()); } else { - gmm::mult_by_row_bwd(matrix, x, x, gmm::abstract_dense()); + gmm::mult_by_row_bwd(gmmMatrix, x, x, gmm::abstract_dense()); } } - template - void GmmxxMultiplier::multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { - std::vector* target = &result; - std::unique_ptr> temporary; + template + void GmmxxMultiplier::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) const { + initialize(); + std::vector* target = &result; if (&x == &result) { - STORM_LOG_WARN("Using temporary in 'multAddReduce'."); - temporary = std::make_unique>(x.size()); - target = temporary.get(); + if (this->cachedVector) { + this->cachedVector->resize(x.size()); + } else { + this->cachedVector = std::make_unique>(x.size()); + } + target = this->cachedVector.get(); } - - if (this->parallelize()) { - multAddReduceParallel(dir, rowGroupIndices, matrix, x, b, *target, choices); + if (parallelize(env)) { + multAddReduceParallel(dir, rowGroupIndices, x, b, *target, choices); } else { - multAddReduceHelper(dir, rowGroupIndices, matrix, x, b, *target, choices); + multAddReduceHelper(dir, rowGroupIndices, x, b, *target, choices); + } + if (&x == &result) { + std::swap(result, *this->cachedVector); } } - template - void GmmxxMultiplier::multAddReduceGaussSeidel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector& x, std::vector const* b, std::vector* choices) const { - multAddReduceHelper(dir, rowGroupIndices, matrix, x, b, x, choices); + 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 { + initialize(); + multAddReduceHelper(dir, rowGroupIndices, x, b, x, choices); } - template - void GmmxxMultiplier::multAddReduceHelper(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { - typedef std::vector VectorType; - typedef gmm::csr_matrix MatrixType; + template + ValueType GmmxxMultiplier::multiplyRow(uint64_t const& rowIndex, std::vector const& x, ValueType const& offset) const { + initialize(); + return vect_sp(gmm::mat_const_row(gmmMatrix, rowIndex), x, typename gmm::linalg_traits>::storage_type(), typename gmm::linalg_traits>::storage_type()) + offset; + } + + template + void GmmxxMultiplier::multAdd(std::vector const& x, std::vector const* b, std::vector& result) const { + if (b) { + gmm::mult_add(gmmMatrix, x, *b, result); + } else { + gmm::mult(gmmMatrix, x, result); + } + } + + template + void GmmxxMultiplier::multAddReduceHelper(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { + typedef std::vector VectorType; + typedef gmm::csr_matrix MatrixType; typename gmm::linalg_traits::const_iterator add_it, add_ite; if (b) { @@ -72,7 +128,7 @@ namespace storm { add_ite = gmm::vect_begin(*b) - 1; } typename gmm::linalg_traits::iterator target_it = gmm::vect_end(result) - 1; - typename gmm::linalg_traits::const_row_iterator itr = mat_row_const_end(matrix) - 1; + typename gmm::linalg_traits::const_row_iterator itr = mat_row_const_end(gmmMatrix) - 1; typename std::vector::iterator choice_it; if (choices) { choice_it = choices->end() - 1; @@ -84,7 +140,7 @@ namespace storm { *choice_it = 0; } - T currentValue = storm::utility::zero(); + ValueType currentValue = storm::utility::zero(); // Only multiply and reduce if the row group is not empty. if (*row_group_it != *(row_group_it + 1)) { @@ -103,7 +159,7 @@ namespace storm { --itr; for (uint64_t row = *row_group_it + 1, rowEnd = *(row_group_it + 1); row < rowEnd; ++row, --itr, --add_it) { - T newValue = b ? *add_it : storm::utility::zero(); + ValueType newValue = b ? *add_it : storm::utility::zero(); newValue += vect_sp(gmm::linalg_traits::row(itr), x); if (choices) { @@ -127,29 +183,29 @@ namespace storm { } template<> - void GmmxxMultiplier::multAddReduceHelper(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { + void GmmxxMultiplier::multAddReduceHelper(storm::solver::OptimizationDirection const& dir, 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."); } - template - void GmmxxMultiplier::multAddParallel(gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result) const { + template + void GmmxxMultiplier::multAddParallel(std::vector const& x, std::vector const* b, std::vector& result) const { #ifdef STORM_HAVE_INTELTBB if (b) { - gmm::mult_add_parallel(matrix, x, *b, result); + gmm::mult_add_parallel(x, *b, result); } else { - gmm::mult_parallel(matrix, x, result); + gmm::mult_parallel(x, result); } #else STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version."); - multAdd(matrix, x, b, result); + multAdd(x, b, result); #endif } #ifdef STORM_HAVE_INTELTBB - template + template class TbbMultAddReduceFunctor { public: - TbbMultAddReduceFunctor(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) : dir(dir), rowGroupIndices(rowGroupIndices), matrix(matrix), x(x), b(b), result(result), choices(choices) { + TbbMultAddReduceFunctor(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) : dir(dir), rowGroupIndices(rowGroupIndices), matrix(matrix), x(x), b(b), result(result), choices(choices) { // Intentionally left empty. } @@ -158,7 +214,7 @@ namespace storm { auto groupIte = rowGroupIndices.begin() + range.end(); auto itr = mat_row_const_begin(matrix) + *groupIt; - typename std::vector::const_iterator bIt; + typename std::vector::const_iterator bIt; if (b) { bIt = b->begin() + *groupIt; } @@ -174,7 +230,7 @@ namespace storm { *choiceIt = 0; } - T currentValue = storm::utility::zero(); + ValueType currentValue = storm::utility::zero(); // Only multiply and reduce if the row group is not empty. if (*groupIt != *(groupIt + 1)) { @@ -186,7 +242,7 @@ namespace storm { ++itr; for (auto itre = mat_row_const_begin(matrix) + *(groupIt + 1); itr != itre; ++itr) { - T newValue = vect_sp(gmm::linalg_traits>::row(itr), x, typename gmm::linalg_traits>::storage_type(), typename gmm::linalg_traits>::storage_type()); + ValueType newValue = vect_sp(gmm::linalg_traits>::row(itr), x, typename gmm::linalg_traits>::storage_type(), typename gmm::linalg_traits>::storage_type()); if (b) { newValue += *bIt; ++bIt; @@ -208,37 +264,29 @@ namespace storm { private: storm::solver::OptimizationDirection dir; std::vector const& rowGroupIndices; - gmm::csr_matrix const& matrix; - std::vector const& x; - std::vector const* b; - std::vector& result; + gmm::csr_matrix const& matrix; + std::vector const& x; + std::vector const* b; + std::vector& result; std::vector* choices; }; #endif - template - void GmmxxMultiplier::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { + template + void GmmxxMultiplier::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { #ifdef STORM_HAVE_INTELTBB - tbb::parallel_for(tbb::blocked_range(0, rowGroupIndices.size() - 1, 10), TbbMultAddReduceFunctor(dir, rowGroupIndices, matrix, x, b, result, choices)); + tbb::parallel_for(tbb::blocked_range(0, rowGroupIndices.size() - 1, 10), TbbMultAddReduceFunctor(dir, rowGroupIndices, this->gmmMatrix, x, b, result, choices)); #else STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version."); - multAddReduce(dir, rowGroupIndices, matrix, x, b, result, choices); + multAddReduceHelper(dir, rowGroupIndices, x, b, result, choices); #endif } template<> - void GmmxxMultiplier::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { + void GmmxxMultiplier::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, 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, "This operation is not supported."); } - template - T GmmxxMultiplier::multiplyRow(gmm::csr_matrix const& matrix, uint64_t const& rowIndex, std::vector const& x) const { - return vect_sp(gmm::mat_const_row(matrix, rowIndex), x, typename gmm::linalg_traits>::storage_type(), typename gmm::linalg_traits>::storage_type()); - } - - - - template class GmmxxMultiplier; #ifdef STORM_HAVE_CARL diff --git a/src/storm/solver/GmmxxMultiplier.h b/src/storm/solver/GmmxxMultiplier.h index aebdb069d..9ea4089b9 100644 --- a/src/storm/solver/GmmxxMultiplier.h +++ b/src/storm/solver/GmmxxMultiplier.h @@ -1,32 +1,42 @@ #pragma once -#include "storm/utility/VectorHelper.h" +#include "storm/solver/Multiplier.h" #include "storm/adapters/GmmxxAdapter.h" #include "storm-config.h" namespace storm { + + namespace storage { + template + class SparseMatrix; + } + namespace solver { - template - class GmmxxMultiplier : public storm::utility::VectorHelper { + template + class GmmxxMultiplier : public Multiplier { public: - GmmxxMultiplier(); + GmmxxMultiplier(storm::storage::SparseMatrix const& matrix); - void multAdd(gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result) const; - void multAddGaussSeidelBackward(gmm::csr_matrix const& matrix, std::vector& x, std::vector const* b) const; + 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 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 ValueType multiplyRow(uint64_t const& rowIndex, std::vector const& x, ValueType const& offset) const override; + virtual void clearCache() const override; + private: + void initialize() const; - void multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const; - void multAddReduceGaussSeidel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector& x, std::vector const* b, std::vector* choices = nullptr) const; + bool parallelize(Environment const& env) const; - void multAddParallel(gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result) const; - void multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const; - - T multiplyRow(gmm::csr_matrix const& matrix, uint64_t const& rowIndex, std::vector const& x) const; - - private: - void multAddReduceHelper(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const; + 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; + + mutable gmm::csr_matrix gmmMatrix; }; }