diff --git a/resources/3rdparty/gmm-5.2/include/gmm/gmm_blas.h b/resources/3rdparty/gmm-5.2/include/gmm/gmm_blas.h index de19b243d..19e58f755 100644 --- a/resources/3rdparty/gmm-5.2/include/gmm/gmm_blas.h +++ b/resources/3rdparty/gmm-5.2/include/gmm/gmm_blas.h @@ -1706,7 +1706,7 @@ namespace gmm { {} }; - template + template class tbbHelper_mult_add_by_row { L2 const* my_l2; L3 const* my_l3; @@ -1714,6 +1714,7 @@ namespace gmm { // Typedefs for Iterator Types typedef typename linalg_traits::iterator frm_IT1; typedef typename linalg_traits::const_row_iterator frm_IT2; + typedef typename linalg_traits::const_iterator frm_IT3; public: void operator()( const forward_range_mult& r ) const { @@ -1722,7 +1723,7 @@ namespace gmm { frm_IT1 it = r.begin(); frm_IT1 ite = r.end(); frm_IT2 itr = r.begin_row(); - frm_IT1 addIt = my_l3->begin(); + frm_IT3 addIt = my_l3->begin(); for (; it != ite; ++it, ++itr, ++addIt) { *it = vect_sp(linalg_traits::row(itr), l2, @@ -1736,14 +1737,65 @@ namespace gmm { } }; + + template + class TbbMultFunctor { + public: + TbbMultFunctor(L1 const& l1, L2 const& l2, L3& l3) : l1(l1), l2(l2), l3(l3) { + // Intentionally left empty. + } + + void operator()(tbb::blocked_range const& range) const { + auto itr = mat_row_const_begin(l1) + range.begin(); + auto l2it = l2.begin() + range.begin(); + auto l3it = l3.begin() + range.begin(); + auto l3ite = l3.begin() + range.end(); + + for (; l3it != l3ite; ++l3it, ++l2it, ++itr) { + *l3it = vect_sp(linalg_traits::row(itr), l2, typename linalg_traits::storage_type(), typename linalg_traits::storage_type()); + } + } + + private: + L1 const& l1; + L2 const& l2; + L3& l3; + }; + template void mult_by_row_parallel(const L1& l1, const L2& l2, L3& l3, abstract_dense) { - tbb::parallel_for(forward_range_mult::iterator, typename linalg_traits::const_row_iterator>(it, ite, itr), tbbHelper_mult_by_row(&l2)); + tbb::parallel_for(tbb::blocked_range(0, vect_size(l3), 10), TbbMultFunctor(l1, l2, l3)); } + template + class TbbMultAddFunctor { + public: + TbbMultAddFunctor(L1 const& l1, L2 const& l2, L3 const& l3, L4& l4) : l1(l1), l2(l2), l3(l3), l4(l4) { + // Intentionally left empty. + } + + void operator()(tbb::blocked_range const& range) const { + auto itr = mat_row_const_begin(l1) + range.begin(); + auto l2it = l2.begin() + range.begin(); + auto l3it = l3.begin() + range.begin(); + auto l4it = l4.begin() + range.begin(); + auto l4ite = l4.begin() + range.end(); + + for (; l4it != l4ite; ++l4it, ++l3it, ++l2it, ++itr) { + *l4it = vect_sp(linalg_traits::row(itr), l2, typename linalg_traits::storage_type(), typename linalg_traits::storage_type()) + *l3it; + } + } + + private: + L1 const& l1; + L2 const& l2; + L3 const& l3; + L4& l4; + }; + template void mult_add_by_row_parallel(const L1& l1, const L2& l2, const L3& l3, L4& l4, abstract_dense) { - tbb::parallel_for(forward_range_mult::iterator, typename linalg_traits::const_row_iterator>(it, ite, itr), tbbHelper_mult_add_by_row(&l2, &l3)); + tbb::parallel_for(tbb::blocked_range(0, vect_size(l4), 10), TbbMultAddFunctor(l1, l2, l3, l4)); } #endif @@ -1879,6 +1931,24 @@ namespace gmm { } #ifdef STORM_HAVE_INTELTBB + /** Multiply. l3 = l1*l2; */ + template inline + void mult_parallel(const L1& l1, const L2& l2, L3& l3) { + size_type m = mat_nrows(l1), n = mat_ncols(l1); + if (!m || !n) return; + GMM_ASSERT2(n==vect_size(l2), "dimensions mismatch"); + if (!same_origin(l2, l3)) { + mult_parallel_spec(l1, l2, l3, typename principal_orientation_type::sub_orientation>::potype()); + } else { + GMM_WARNING2("Warning, temporaries are used for mult\n"); + typename temporary_vector::vector_type l2tmp(vect_size(l2)); + copy(l2, l2tmp); + mult_parallel_spec(l1, l2tmp, l3, typename principal_orientation_type::sub_orientation>::potype()); + } + } + /** Multiply-accumulate. l4 = l1*l2 + l3; */ template inline void mult_add_parallel(const L1& l1, const L2& l2, const L3& l3, L4& l4) { @@ -1886,16 +1956,14 @@ namespace gmm { if (!m || !n) return; GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l3) && vect_size(l3) == vect_size(l4), "dimensions mismatch"); if (!same_origin(l2, l4) && !same_origin(l3, l4) && !same_origin(l2, l3)) { - mult_add_parallel_spec(l1, l2, l3, l4, typename principal_orientation_type::sub_orientation>::potype()); + mult_add_parallel_spec(l1, l2, l3, l4, typename principal_orientation_type::sub_orientation>::potype()); } else { GMM_WARNING2("Warning, Multiple temporaries are used for mult\n"); typename temporary_vector::vector_type l2tmp(vect_size(l2)); copy(l2, l2tmp); typename temporary_vector::vector_type l3tmp(vect_size(l3)); copy(l3, l3tmp); - mult_add_parallel_spec(l1, l2tmp, l3tmp, l4, typename principal_orientation_type::sub_orientation>::potype()); + mult_add_parallel_spec(l1, l2tmp, l3tmp, l4, typename principal_orientation_type::sub_orientation>::potype()); } } #endif diff --git a/src/storm/adapters/GmmxxAdapter.cpp b/src/storm/adapters/GmmxxAdapter.cpp index 63769cd34..4210a00ee 100644 --- a/src/storm/adapters/GmmxxAdapter.cpp +++ b/src/storm/adapters/GmmxxAdapter.cpp @@ -5,8 +5,7 @@ #include "storm/adapters/RationalNumberAdapter.h" #include "storm/adapters/RationalFunctionAdapter.h" -#include "storm/utility/constants.h" -#include "storm/exceptions/NotSupportedException.h" +#include "storm/utility/macros.h" namespace storm { namespace adapters { @@ -42,125 +41,12 @@ namespace storm { return result; } - - template - void GmmxxMultiplier::multAdd(gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result) { - if (b) { - gmm::mult_add(matrix, x, *b, result); - } else { - gmm::mult(matrix, x, result); - } - } - - template - void GmmxxMultiplier::multAddGaussSeidelBackward(gmm::csr_matrix const& matrix, std::vector& x, std::vector const* b) { - STORM_LOG_ASSERT(matrix.nr == matrix.nc, "Expecting square matrix."); - if (b) { - gmm::mult_add_by_row_bwd(matrix, x, *b, x, gmm::abstract_dense()); - } else { - gmm::mult_by_row_bwd(matrix, 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) { - std::vector* target = &result; - std::unique_ptr> temporary; - if (&x == &result) { - STORM_LOG_WARN("Using temporary in 'multAddReduce'."); - temporary = std::make_unique>(x.size()); - target = temporary.get(); - } - - multAddReduceHelper(dir, rowGroupIndices, matrix, x, b, *target, choices); - } - - 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) { - multAddReduceHelper(dir, rowGroupIndices, matrix, 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) { - typedef std::vector VectorType; - typedef gmm::csr_matrix MatrixType; - - 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; - } - 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 std::vector::iterator choice_it; - if (choices) { - choice_it = choices->end() - 1; - } - - uint64_t choice; - 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) { - T currentValue = b ? *add_it : storm::utility::zero(); - currentValue += vect_sp(gmm::linalg_traits::row(itr), x); - - if (choices) { - choice = *(row_group_it + 1) - 1 - *row_group_it; - *choice_it = choice; - } - - --itr; - if (b) { - --add_it; - } - - for (uint64_t row = *row_group_it + 1, rowEnd = *(row_group_it + 1); row < rowEnd; ++row, --itr) { - T newValue = b ? *add_it : storm::utility::zero(); - newValue += vect_sp(gmm::linalg_traits::row(itr), x); - - if (choices) { - --choice; - } - - if ((dir == OptimizationDirection::Minimize && newValue < currentValue) || (dir == OptimizationDirection::Maximize && newValue > currentValue)) { - currentValue = newValue; - if (choices) { - *choice_it = choice; - } - } - if (b) { - --add_it; - } - } - - // Write back final value. - *target_it = currentValue; - } - } - - 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) { - STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Operation not supported for this data type."); - } - -#ifdef STORM_HAVE_INTELTBB - template - void GmmxxMultiplier::multAddParallel(gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result) { - if (b) { - gmm::mult_add_parallel(matrix, x, *b, result); - } else { - gmm::mult_parallel(matrix, x, result); - } - } -#endif template class GmmxxAdapter; - template class GmmxxMultiplier; #ifdef STORM_HAVE_CARL template class GmmxxAdapter; template class GmmxxAdapter; - - template class GmmxxMultiplier; - template class GmmxxMultiplier; #endif } diff --git a/src/storm/adapters/GmmxxAdapter.h b/src/storm/adapters/GmmxxAdapter.h index 442b2575f..a5e139db0 100644 --- a/src/storm/adapters/GmmxxAdapter.h +++ b/src/storm/adapters/GmmxxAdapter.h @@ -6,8 +6,6 @@ #include "storm/storage/SparseMatrix.h" -#include "storm/utility/macros.h" - namespace storm { namespace adapters { @@ -21,22 +19,5 @@ namespace storm { static std::unique_ptr> toGmmxxSparseMatrix(storm::storage::SparseMatrix const& matrix); }; - template - class GmmxxMultiplier { - public: - static void multAdd(gmm::csr_matrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result); - static void multAddGaussSeidelBackward(gmm::csr_matrix const& matrix, std::vector& x, std::vector const* b); - - static 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); - static 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); - -#ifdef STORM_HAVE_INTELTBB - static void multAddParallel(gmm::csr_matrix const& matrix, std::vector const& x, std::vector const& b, std::vector& result); -#endif - - private: - static 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); - }; - } } diff --git a/src/storm/settings/modules/CoreSettings.cpp b/src/storm/settings/modules/CoreSettings.cpp index 234a1b9a1..836943dae 100644 --- a/src/storm/settings/modules/CoreSettings.cpp +++ b/src/storm/settings/modules/CoreSettings.cpp @@ -12,6 +12,7 @@ #include "storm/utility/macros.h" #include "storm/exceptions/IllegalArgumentValueException.h" +#include "storm/exceptions/InvalidOptionException.h" namespace storm { namespace settings { @@ -31,6 +32,8 @@ namespace storm { const std::string CoreSettings::engineOptionShortName = "e"; const std::string CoreSettings::ddLibraryOptionName = "ddlib"; const std::string CoreSettings::cudaOptionName = "cuda"; + const std::string CoreSettings::intelTbbOptionName = "enable-tbb"; + const std::string CoreSettings::intelTbbOptionShortName = "tbb"; CoreSettings::CoreSettings() : ModuleSettings(moduleName), engine(CoreSettings::Engine::Sparse) { this->addOption(storm::settings::OptionBuilder(moduleName, counterexampleOptionName, false, "Generates a counterexample for the given PRCTL formulas if not satisfied by the model.").setShortName(counterexampleOptionShortName).build()); @@ -56,7 +59,9 @@ namespace storm { this->addOption(storm::settings::OptionBuilder(moduleName, smtSolverOptionName, false, "Sets which SMT solver is preferred.") .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of an SMT solver.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(smtSolvers)).setDefaultValueString("z3").build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, statisticsOptionName, false, "Sets whether to display statistics if available.").setShortName(statisticsOptionShortName).build()); + this->addOption(storm::settings::OptionBuilder(moduleName, cudaOptionName, false, "Sets whether to use CUDA.").build()); + this->addOption(storm::settings::OptionBuilder(moduleName, intelTbbOptionName, false, "Sets whether to use Intel TBB (if Storm was built with support for TBB).").setShortName(intelTbbOptionShortName).build()); } bool CoreSettings::isCounterexampleSet() const { @@ -127,8 +132,12 @@ namespace storm { bool CoreSettings::isShowStatisticsSet() const { return this->getOption(statisticsOptionName).getHasOptionBeenSet(); } - - bool CoreSettings::isCudaSet() const { + + bool CoreSettings::isUseIntelTbbSet() const { + return this->getOption(intelTbbOptionName).getHasOptionBeenSet(); + } + + bool CoreSettings::isUseCudaSet() const { return this->getOption(cudaOptionName).getHasOptionBeenSet(); } @@ -159,7 +168,12 @@ namespace storm { } bool CoreSettings::check() const { +#ifdef STORM_HAVE_INTELTBB + return true; +#else + STORM_LOG_WARN_COND(!isUseIntelTbbSet(), "Enabling TBB is not supported in this version of Storm as it was not built with support for it."); return true; +#endif } } // namespace modules diff --git a/src/storm/settings/modules/CoreSettings.h b/src/storm/settings/modules/CoreSettings.h index b08c01028..1880ee109 100644 --- a/src/storm/settings/modules/CoreSettings.h +++ b/src/storm/settings/modules/CoreSettings.h @@ -116,12 +116,19 @@ namespace storm { */ bool isShowStatisticsSet() const; + /*! + * Retrieves whether the option to use Intel TBB is set. + * + * @return True iff the option was set. + */ + bool isUseIntelTbbSet() const; + /*! * Retrieves whether the option to use CUDA is set. * * @return True iff the option was set. */ - bool isCudaSet() const; + bool isUseCudaSet() const; /*! * Retrieves the selected engine. @@ -157,6 +164,8 @@ namespace storm { static const std::string engineOptionName; static const std::string engineOptionShortName; static const std::string ddLibraryOptionName; + static const std::string intelTbbOptionName; + static const std::string intelTbbOptionShortName; static const std::string cudaOptionName; }; diff --git a/src/storm/solver/GmmxxLinearEquationSolver.cpp b/src/storm/solver/GmmxxLinearEquationSolver.cpp index 26cf7d44a..b922a1b9c 100644 --- a/src/storm/solver/GmmxxLinearEquationSolver.cpp +++ b/src/storm/solver/GmmxxLinearEquationSolver.cpp @@ -4,6 +4,9 @@ #include #include "storm/adapters/GmmxxAdapter.h" + +#include "storm/solver/GmmxxMultiplier.h" + #include "storm/settings/SettingsManager.h" #include "storm/utility/vector.h" #include "storm/utility/constants.h" @@ -187,7 +190,7 @@ namespace storm { template void GmmxxLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { - storm::adapters::GmmxxMultiplier::multAdd(*gmmxxA, x, b, result); + multiplier.multAdd(*gmmxxA, x, b, result); if (!this->isCachingEnabled()) { clearCache(); @@ -196,7 +199,7 @@ namespace storm { template void GmmxxLinearEquationSolver::multiplyAndReduce(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector& result, std::vector* choices) const { - storm::adapters::GmmxxMultiplier::multAddReduce(dir, rowGroupIndices, *gmmxxA, x, b, result, choices); + multiplier.multAddReduce(dir, rowGroupIndices, *gmmxxA, x, b, result, choices); } template @@ -206,12 +209,12 @@ namespace storm { template void GmmxxLinearEquationSolver::multiplyGaussSeidel(std::vector& x, std::vector const* b) const { - storm::adapters::GmmxxMultiplier::multAddGaussSeidelBackward(*gmmxxA, x, b); + multiplier.multAddGaussSeidelBackward(*gmmxxA, x, b); } template void GmmxxLinearEquationSolver::multiplyAndReduceGaussSeidel(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices) const { - storm::adapters::GmmxxMultiplier::multAddReduceGaussSeidel(dir, rowGroupIndices, *gmmxxA, x, b, choices); + multiplier.multAddReduceGaussSeidel(dir, rowGroupIndices, *gmmxxA, x, b, choices); } template diff --git a/src/storm/solver/GmmxxLinearEquationSolver.h b/src/storm/solver/GmmxxLinearEquationSolver.h index fe301de12..cce3a6369 100644 --- a/src/storm/solver/GmmxxLinearEquationSolver.h +++ b/src/storm/solver/GmmxxLinearEquationSolver.h @@ -5,7 +5,9 @@ #include "storm/utility/gmm.h" -#include "LinearEquationSolver.h" +#include "storm/solver/GmmxxMultiplier.h" + +#include "storm/solver/LinearEquationSolver.h" namespace storm { namespace solver { @@ -95,7 +97,6 @@ namespace storm { void setSettings(GmmxxLinearEquationSolverSettings const& newSettings); GmmxxLinearEquationSolverSettings const& getSettings() const; - virtual void clearCache() const override; private: @@ -108,6 +109,9 @@ namespace storm { // The settings used by the solver. GmmxxLinearEquationSolverSettings settings; + // A multiplier object used to dispatch the multiplication calls. + GmmxxMultiplier multiplier; + // cached data obtained during solving mutable std::unique_ptr>> iluPreconditioner; mutable std::unique_ptr>> diagonalPreconditioner; diff --git a/src/storm/solver/GmmxxMultiplier.cpp b/src/storm/solver/GmmxxMultiplier.cpp new file mode 100644 index 000000000..ec196af49 --- /dev/null +++ b/src/storm/solver/GmmxxMultiplier.cpp @@ -0,0 +1,221 @@ +#include "storm/solver/GmmxxMultiplier.h" + +#include "storm/adapters/RationalNumberAdapter.h" +#include "storm/adapters/RationalFunctionAdapter.h" + +#include "storm/utility/constants.h" +#include "storm/exceptions/NotSupportedException.h" + +#include "storm/utility/macros.h" + +namespace storm { + namespace solver { + + template + GmmxxMultiplier::GmmxxMultiplier() : storm::utility::VectorHelper() { + // 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); + } else { + gmm::mult(matrix, x, result); + } + } + } + + 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."); + if (b) { + gmm::mult_add_by_row_bwd(matrix, x, *b, x, gmm::abstract_dense()); + } else { + gmm::mult_by_row_bwd(matrix, 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; + if (&x == &result) { + STORM_LOG_WARN("Using temporary in 'multAddReduce'."); + temporary = std::make_unique>(x.size()); + target = temporary.get(); + } + + if (this->parallelize()) { + multAddReduceParallel(dir, rowGroupIndices, matrix, x, b, *target, choices); + } else { + multAddReduceHelper(dir, rowGroupIndices, matrix, x, b, *target, choices); + } + } + + 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::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; + + 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; + } + 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 std::vector::iterator choice_it; + if (choices) { + choice_it = choices->end() - 1; + } + + uint64_t choice; + 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) { + T currentValue = b ? *add_it : storm::utility::zero(); + currentValue += vect_sp(gmm::linalg_traits::row(itr), x); + + if (choices) { + choice = *(row_group_it + 1) - 1 - *row_group_it; + *choice_it = choice; + } + + --itr; + if (b) { + --add_it; + } + + for (uint64_t row = *row_group_it + 1, rowEnd = *(row_group_it + 1); row < rowEnd; ++row, --itr) { + T newValue = b ? *add_it : storm::utility::zero(); + newValue += vect_sp(gmm::linalg_traits::row(itr), x); + + if (choices) { + --choice; + } + + if ((dir == OptimizationDirection::Minimize && newValue < currentValue) || (dir == OptimizationDirection::Maximize && newValue > currentValue)) { + currentValue = newValue; + if (choices) { + *choice_it = choice; + } + } + if (b) { + --add_it; + } + } + + // Write back final value. + *target_it = currentValue; + } + } + + 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 { + 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 { +#ifdef STORM_HAVE_INTELTBB + if (b) { + gmm::mult_add_parallel(matrix, x, *b, result); + } else { + gmm::mult_parallel(matrix, x, result); + } +#else + STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version."); + multAdd(matrix, x, b, result); +#endif + } + + 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) { + // Intentionally left empty. + } + + void operator()(tbb::blocked_range const& range) const { + auto groupIt = rowGroupIndices.begin() + range.begin(); + auto groupIte = rowGroupIndices.begin() + range.end(); + + auto itr = mat_row_const_begin(matrix) + *groupIt; + typename std::vector::const_iterator bIt; + if (b) { + bIt = b->begin() + *groupIt; + } + typename std::vector::iterator choiceIt; + if (choices) { + choiceIt = choices->begin() + range.begin(); + } + + auto resultIt = result.begin() + range.begin(); + + for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) { + T currentValue = vect_sp(gmm::linalg_traits>::row(itr), x, typename gmm::linalg_traits>::storage_type(), typename gmm::linalg_traits>::storage_type()); + if (b) { + currentValue += *bIt; + ++bIt; + } + if (choices) { + *choiceIt = 0; + } + + ++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()); + if (b) { + newValue += *bIt; + ++bIt; + } + + if ((dir == OptimizationDirection::Minimize && newValue < currentValue) || (dir == OptimizationDirection::Maximize && newValue > currentValue)) { + currentValue = newValue; + if (choices) { + *choiceIt = std::distance(mat_row_const_begin(matrix), itr) - *groupIt; + } + } + } + + *resultIt = currentValue; + } + } + + 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; + std::vector* choices; + }; + + 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 { + tbb::parallel_for(tbb::blocked_range(0, rowGroupIndices.size() - 1, 10), TbbMultAddReduceFunctor(dir, rowGroupIndices, matrix, x, b, result, choices)); + } + + 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 { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); + } + + template class GmmxxMultiplier; + +#ifdef STORM_HAVE_CARL + template class GmmxxMultiplier; + template class GmmxxMultiplier; +#endif + + } +} diff --git a/src/storm/solver/GmmxxMultiplier.h b/src/storm/solver/GmmxxMultiplier.h new file mode 100644 index 000000000..237f56b63 --- /dev/null +++ b/src/storm/solver/GmmxxMultiplier.h @@ -0,0 +1,31 @@ +#pragma once + +#include "storm/utility/VectorHelper.h" + +#include "storm/adapters/GmmxxAdapter.h" + +#include "storm-config.h" + +namespace storm { + namespace solver { + + template + class GmmxxMultiplier : public storm::utility::VectorHelper { + public: + GmmxxMultiplier(); + + 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; + + 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; + + 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; + + 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; + }; + + } +} diff --git a/src/storm/solver/LinearEquationSolver.cpp b/src/storm/solver/LinearEquationSolver.cpp index 2f7d60287..bdcc4041e 100644 --- a/src/storm/solver/LinearEquationSolver.cpp +++ b/src/storm/solver/LinearEquationSolver.cpp @@ -69,7 +69,7 @@ namespace storm { setCachingEnabled(true); this->multiply(x, b, *cachedRowVector); - storm::utility::vector::reduceVectorMinOrMax(dir, *cachedRowVector, result, rowGroupIndices, choices); + vectorHelper.reduceVector(dir, *cachedRowVector, result, rowGroupIndices, choices); // restore the old caching setting setCachingEnabled(cachingWasEnabled); diff --git a/src/storm/solver/LinearEquationSolver.h b/src/storm/solver/LinearEquationSolver.h index 47e82204a..1d3eb5a83 100644 --- a/src/storm/solver/LinearEquationSolver.h +++ b/src/storm/solver/LinearEquationSolver.h @@ -8,6 +8,8 @@ #include "storm/solver/MultiplicationStyle.h" #include "storm/solver/OptimizationDirection.h" +#include "storm/utility/VectorHelper.h" + #include "storm/storage/SparseMatrix.h" namespace storm { @@ -171,6 +173,9 @@ namespace storm { /// Whether some of the generated data during solver calls should be cached. mutable bool cachingEnabled; + + /// An object that can be used to reduce vectors. + storm::utility::VectorHelper vectorHelper; }; template diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 855b56fed..d1d988e50 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -171,10 +171,11 @@ namespace storm { // Set up additional environment variables. uint_fast64_t iterationCount = 0; bool converged = false; - + while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) { // Compute D^-1 * (b - LU * x) and store result in nextX. - jacobiLU.multiplyWithVector(*currentX, *nextX); + multiplier.multAdd(jacobiLU, *currentX, nullptr, *nextX); + storm::utility::vector::subtractVectors(b, *nextX, *nextX); storm::utility::vector::multiplyVectorsPointwise(jacobiD, *nextX, *nextX); @@ -292,7 +293,7 @@ namespace storm { // Create a vector that always holds Ax. std::vector currentAx(x.size()); - walkerChaeData->matrix.multiplyWithVector(*currentX, currentAx); + multiplier.multAdd(walkerChaeData->matrix, *currentX, nullptr, currentAx); // (3) Perform iterations until convergence. bool converged = false; @@ -302,7 +303,7 @@ namespace storm { walkerChaeData->matrix.performWalkerChaeStep(*currentX, walkerChaeData->columnSums, walkerChaeData->b, currentAx, *nextX); // Compute new Ax. - walkerChaeData->matrix.multiplyWithVector(*nextX, currentAx); + multiplier.multAdd(walkerChaeData->matrix, *nextX, nullptr, currentAx); // Check for convergence. converged = storm::utility::vector::computeSquaredNorm2Difference(currentAx, walkerChaeData->b) <= squaredErrorBound; @@ -355,14 +356,14 @@ namespace storm { template void NativeLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { if (&x != &result) { - A->multiplyWithVector(x, result, b); + multiplier.multAdd(*A, x, b, result); } else { // If the two vectors are aliases, we need to create a temporary. if (!this->cachedRowVector) { this->cachedRowVector = std::make_unique>(getMatrixRowCount()); } - A->multiplyWithVector(x, *this->cachedRowVector, b); + multiplier.multAdd(*A, x, b, *this->cachedRowVector); result.swap(*this->cachedRowVector); if (!this->isCachingEnabled()) { @@ -374,14 +375,14 @@ namespace storm { template void NativeLinearEquationSolver::multiplyAndReduce(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector& result, std::vector* choices) const { if (&x != &result) { - A->multiplyAndReduce(dir, rowGroupIndices, x, b, result, choices); + multiplier.multAddReduce(dir, rowGroupIndices, *A, x, b, result, choices); } else { // If the two vectors are aliases, we need to create a temporary. if (!this->cachedRowVector) { this->cachedRowVector = std::make_unique>(getMatrixRowCount()); } - this->A->multiplyAndReduce(dir, rowGroupIndices, x, b, *this->cachedRowVector, choices); + multiplier.multAddReduce(dir, rowGroupIndices, *A, x, b, *this->cachedRowVector, choices); result.swap(*this->cachedRowVector); if (!this->isCachingEnabled()) { @@ -398,12 +399,12 @@ namespace storm { template void NativeLinearEquationSolver::multiplyGaussSeidel(std::vector& x, std::vector const* b) const { STORM_LOG_ASSERT(this->A->getRowCount() == this->A->getColumnCount(), "This function is only applicable for square matrices."); - A->multiplyWithVector(x, x, b, true, storm::storage::SparseMatrix::MultiplicationDirection::Backward); + multiplier.multAddGaussSeidelBackward(*A, x, b); } template void NativeLinearEquationSolver::multiplyAndReduceGaussSeidel(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector& x, std::vector const* b, std::vector* choices) const { - A->multiplyAndReduce(dir, rowGroupIndices, x, b, x, choices, true, storm::storage::SparseMatrix::MultiplicationDirection::Backward); + multiplier.multAddReduceGaussSeidelBackward(dir, rowGroupIndices, *A, x, b, choices); } template diff --git a/src/storm/solver/NativeLinearEquationSolver.h b/src/storm/solver/NativeLinearEquationSolver.h index 2c9633c55..facf02868 100644 --- a/src/storm/solver/NativeLinearEquationSolver.h +++ b/src/storm/solver/NativeLinearEquationSolver.h @@ -3,7 +3,9 @@ #include -#include "LinearEquationSolver.h" +#include "storm/solver/LinearEquationSolver.h" + +#include "storm/solver/NativeMultiplier.h" namespace storm { namespace solver { @@ -76,10 +78,13 @@ namespace storm { // A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix // the pointer refers to localA. storm::storage::SparseMatrix const* A; - + // The settings used by the solver. NativeLinearEquationSolverSettings settings; - + + // An object to dispatch all multiplication operations. + NativeMultiplier multiplier; + // cached auxiliary data mutable std::unique_ptr, std::vector>> jacobiDecomposition; diff --git a/src/storm/solver/NativeMultiplier.cpp b/src/storm/solver/NativeMultiplier.cpp new file mode 100644 index 000000000..21da2ff86 --- /dev/null +++ b/src/storm/solver/NativeMultiplier.cpp @@ -0,0 +1,101 @@ +#include "storm/solver/NativeMultiplier.h" + +#include "storm-config.h" + +#include "storm/storage/SparseMatrix.h" + +#include "storm/adapters/RationalNumberAdapter.h" +#include "storm/adapters/RationalFunctionAdapter.h" + +#include "storm/utility/macros.h" + +namespace storm { + namespace solver { + + template + NativeMultiplier::NativeMultiplier() : storm::utility::VectorHelper() { + // Intentionally left empty. + } + + template + void NativeMultiplier::multAdd(storm::storage::SparseMatrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result) const { + std::vector* target = &result; + std::unique_ptr> temporary; + if (&x == &result) { + STORM_LOG_WARN("Using temporary in 'multAdd'."); + temporary = std::make_unique>(x.size()); + target = temporary.get(); + } + + if (this->parallelize()) { + multAddParallel(matrix, x, b, result); + } else { + matrix.multiplyWithVector(x, result, b); + } + + if (target == temporary.get()) { + std::swap(result, *temporary); + } + } + + template + void NativeMultiplier::multAddGaussSeidelBackward(storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector const* b) const { + matrix.multiplyWithVectorBackward(x, x, b); + } + + template + void NativeMultiplier::multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, storm::storage::SparseMatrix 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; + if (&x == &result) { + STORM_LOG_WARN("Using temporary in 'multAddReduce'."); + temporary = std::make_unique>(x.size()); + target = temporary.get(); + } + + if (this->parallelize()) { + multAddReduceParallel(dir, rowGroupIndices, matrix, x, b, *target, choices); + } else { + matrix.multiplyAndReduce(dir, rowGroupIndices, x, b, *target, choices); + } + + if (target == temporary.get()) { + std::swap(result, *temporary); + } + } + + template + void NativeMultiplier::multAddReduceGaussSeidelBackward(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector const* b, std::vector* choices) const { + matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices); + } + + template + void NativeMultiplier::multAddParallel(storm::storage::SparseMatrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result) const { +#ifdef STORM_HAVE_INTELTBB + matrix.multiplyWithVectorParallel(x, result, b); +#else + STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version."); + multAdd(matrix, x, b, result); +#endif + } + + template + void NativeMultiplier::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, storm::storage::SparseMatrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices) const { +#ifdef STORM_HAVE_INTELTBB + matrix.multiplyAndReduceParallel(dir, rowGroupIndices, x, b, result, choices); +#else + STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version."); + multAddReduce(dir, rowGroupIndices, x, b, result, choices); +#endif + } + + + template class NativeMultiplier; + +#ifdef STORM_HAVE_CARL + template class NativeMultiplier; + template class NativeMultiplier; +#endif + + } +} diff --git a/src/storm/solver/NativeMultiplier.h b/src/storm/solver/NativeMultiplier.h new file mode 100644 index 000000000..01233b1ba --- /dev/null +++ b/src/storm/solver/NativeMultiplier.h @@ -0,0 +1,31 @@ +#pragma once + +#include "storm/utility/VectorHelper.h" + +#include "storm/solver/OptimizationDirection.h" + +namespace storm { + namespace storage { + template + class SparseMatrix; + } + + namespace solver { + + template + class NativeMultiplier : public storm::utility::VectorHelper { + public: + NativeMultiplier(); + + void multAdd(storm::storage::SparseMatrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result) const; + void multAddGaussSeidelBackward(storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector const* b) const; + + void multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, storm::storage::SparseMatrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const; + void multAddReduceGaussSeidelBackward(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector const* b, std::vector* choices = nullptr) const; + + void multAddParallel(storm::storage::SparseMatrix 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, storm::storage::SparseMatrix const& matrix, std::vector const& x, std::vector const* b, std::vector& result, std::vector* choices = nullptr) const; + }; + + } +} diff --git a/src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp b/src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp index 6d70c7ea9..3307e8af4 100644 --- a/src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp @@ -25,7 +25,7 @@ namespace storm { template TopologicalMinMaxLinearEquationSolver::TopologicalMinMaxLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), relative(relative) { // Get the settings object to customize solving. - this->enableCuda = storm::settings::getModule().isCudaSet(); + this->enableCuda = storm::settings::getModule().isUseCudaSet(); #ifdef STORM_HAVE_CUDA STORM_LOG_INFO_COND(this->enableCuda, "Option CUDA was not set, but the topological value iteration solver will use it anyways."); #endif diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index c622ed7d9..258493bd6 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -1304,36 +1304,23 @@ namespace storm { } template - void SparseMatrix::multiplyWithVector(std::vector const& vector, std::vector& result, std::vector const* summand, bool allowAliasing, MultiplicationDirection const& multiplicationDirection) const { + void SparseMatrix::multiplyWithVector(std::vector const& vector, std::vector& result, std::vector const* summand) const { // If the vector and the result are aliases and this is not set to be allowed, we need and temporary vector. std::vector* target; std::vector temporary; - bool vectorsAliased = &vector == &result; - if (!allowAliasing && vectorsAliased) { - STORM_LOG_WARN("Vectors are aliased but are not allowed to be. Using temporary, which is potentially slow."); + if (&vector == &result) { + STORM_LOG_WARN("Vectors are aliased. Using temporary, which is potentially slow."); temporary = std::vector(vector.size()); target = &temporary; - STORM_LOG_WARN_COND(multiplicationDirection != MultiplicationDirection::DontCare, "Not specifying multiplication direction for aliased vectors may yield unexpected results."); } else { target = &result; } - STORM_LOG_WARN_COND(vectorsAliased || multiplicationDirection == MultiplicationDirection::DontCare, "Setting a multiplication direction for unaliased vectors. Check whether this is intended."); - -#ifdef STORM_HAVE_INTELTBB - bool useParallel = !allowAliasing && multiplicationDirection == MultiplicationDirection::DontCare && this->getNonzeroEntryCount() > 10000; - if (useParallel) { - this->multiplyWithVectorParallel(vector, *target, summand); - } else { -#endif - if (multiplicationDirection == MultiplicationDirection::Forward || (multiplicationDirection == MultiplicationDirection::DontCare && !vectorsAliased)) { - this->multiplyWithVectorForward(vector, *target, summand); - } else { - this->multiplyWithVectorBackward(vector, *target, summand); - } -#ifdef STORM_HAVE_INTELTBB + this->multiplyWithVectorForward(vector, *target, summand); + + if (target == &temporary) { + std::swap(result, *target); } -#endif } template @@ -1389,6 +1376,47 @@ namespace storm { } #ifdef STORM_HAVE_INTELTBB + template + class TbbMultAddFunctor { + public: + typedef typename storm::storage::SparseMatrix::index_type index_type; + typedef typename storm::storage::SparseMatrix::value_type value_type; + typedef typename storm::storage::SparseMatrix::const_iterator const_iterator; + + TbbMultAddFunctor(std::vector> const& columnsAndEntries, std::vector const& rowIndications, std::vector const& x, std::vector& result, std::vector const* summand) : columnsAndEntries(columnsAndEntries), rowIndications(rowIndications), x(x), result(result), summand(summand) { + // Intentionally left empty. + } + + void operator()(tbb::blocked_range const& range) const { + index_type startRow = range.begin(); + index_type endRow = range.end(); + typename std::vector::const_iterator rowIterator = rowIndications.begin() + startRow; + const_iterator it = columnsAndEntries.begin() + *rowIterator; + const_iterator ite; + typename std::vector::iterator resultIterator = result.begin() + startRow; + typename std::vector::iterator resultIteratorEnd = result.begin() + endRow; + typename std::vector::const_iterator summandIterator; + if (summand) { + summandIterator = summand->begin() + startRow; + } + + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++summandIterator) { + *resultIterator = summand ? *summandIterator : storm::utility::zero(); + + for (ite = columnsAndEntries.begin() + *(rowIterator + 1); it != ite; ++it) { + *resultIterator += it->getValue() * x[it->getColumn()]; + } + } + } + + private: + std::vector> const& columnsAndEntries; + std::vector const& rowIndications; + std::vector const& x; + std::vector& result; + std::vector const* summand; + }; + template void SparseMatrix::multiplyWithVectorParallel(std::vector const& vector, std::vector& result, std::vector const* summand) const { if (&vector == &result) { @@ -1397,34 +1425,7 @@ namespace storm { multiplyWithVectorParallel(vector, tmpVector); result = std::move(tmpVector); } else { - tbb::parallel_for(tbb::blocked_range(0, result.size(), 10), - [&] (tbb::blocked_range const& range) { - index_type startRow = range.begin(); - index_type endRow = range.end(); - const_iterator it = this->begin(startRow); - const_iterator ite; - std::vector::const_iterator rowIterator = this->rowIndications.begin() + startRow; - std::vector::const_iterator rowIteratorEnd = this->rowIndications.begin() + endRow; - typename std::vector::iterator resultIterator = result.begin() + startRow; - typename std::vector::iterator resultIteratorEnd = result.begin() + endRow; - typename std::vector::const_iterator summandIterator; - if (summand) { - summandIterator = summand->begin() + startRow; - } - - for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) { - if (summand) { - *resultIterator = *summandIterator; - ++summandIterator; - } else { - *resultIterator = storm::utility::zero(); - } - - for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { - *resultIterator += it->getValue() * vector[it->getColumn()]; - } - } - }); + tbb::parallel_for(tbb::blocked_range(0, result.size(), 10), TbbMultAddFunctor(columnsAndValues, rowIndications, vector, result, summand)); } } #endif @@ -1560,7 +1561,7 @@ namespace storm { #ifdef STORM_HAVE_CARL template<> void SparseMatrix::multiplyAndReduceForward(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const { - STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This operation is not supported."); + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); } #endif @@ -1618,96 +1619,111 @@ namespace storm { #ifdef STORM_HAVE_CARL template<> void SparseMatrix::multiplyAndReduceBackward(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const { - STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This operation is not supported."); + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); } #endif #ifdef STORM_HAVE_INTELTBB + template + class TbbMultAddReduceFunctor { + public: + typedef typename storm::storage::SparseMatrix::index_type index_type; + typedef typename storm::storage::SparseMatrix::value_type value_type; + typedef typename storm::storage::SparseMatrix::const_iterator const_iterator; + + TbbMultAddReduceFunctor(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector> const& columnsAndEntries, std::vector const& rowIndications, std::vector const& x, std::vector& result, std::vector const* summand, std::vector* choices) : dir(dir), rowGroupIndices(rowGroupIndices), columnsAndEntries(columnsAndEntries), rowIndications(rowIndications), x(x), result(result), summand(summand), choices(choices) { + // Intentionally left empty. + } + + void operator()(tbb::blocked_range const& range) const { + auto groupIt = rowGroupIndices.begin() + range.begin(); + auto groupIte = rowGroupIndices.begin() + range.end(); + + auto rowIt = rowIndications.begin() + *groupIt; + auto elementIt = columnsAndEntries.begin() + *rowIt; + typename std::vector::const_iterator summandIt; + if (summand) { + summandIt = summand->begin() + *groupIt; + } + typename std::vector::iterator choiceIt; + if (choices) { + choiceIt = choices->begin() + range.begin(); + } + + auto resultIt = result.begin() + range.begin(); + + for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) { + ValueType currentValue = summand ? *summandIt : storm::utility::zero(); + + for (auto elementIte = columnsAndEntries.begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { + currentValue += elementIt->getValue() * x[elementIt->getColumn()]; + } + if (choices) { + *choiceIt = 0; + } + + ++rowIt; + ++summandIt; + + for (; static_cast(std::distance(rowIndications.begin(), rowIt)) < *(groupIt + 1); ++rowIt, ++summandIt) { + ValueType newValue = summand ? *summandIt : storm::utility::zero(); + for (auto elementIte = columnsAndEntries.begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { + newValue += elementIt->getValue() * x[elementIt->getColumn()]; + } + + if ((dir == OptimizationDirection::Minimize && newValue < currentValue) || (dir == OptimizationDirection::Maximize && newValue > currentValue)) { + currentValue = newValue; + if (choices) { + *choiceIt = std::distance(rowIndications.begin(), rowIt) - *groupIt; + } + } + } + + // Finally write value to target vector. + *resultIt = currentValue; + } + } + + private: + OptimizationDirection dir; + std::vector const& rowGroupIndices; + std::vector> const& columnsAndEntries; + std::vector const& rowIndications; + std::vector const& x; + std::vector& result; + std::vector const* summand; + std::vector* choices; + }; + template void SparseMatrix::multiplyAndReduceParallel(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { - tbb::parallel_for(tbb::blocked_range(0, rowGroupIndices.size() - 1, 10), - [&] (tbb::blocked_range const& range) { - index_type startRowGroup = range.begin(); - index_type endRowGroup = range.end(); - - auto rowGroupIt = rowGroupIndices.begin() + startRowGroup; - auto rowIt = rowIndications.begin() + startRowGroup; - auto elementIt = this->begin(*rowIt); - typename std::vector::const_iterator summandIt; - if (summand) { - summandIt = summand->begin(); - } - typename std::vector::iterator choiceIt; - if (choices) { - choiceIt = choices->begin() + startRowGroup; - } - - for (auto resultIt = result.begin() + startRowGroup, resultIte = result.begin() + endRow; resultIt != resultIte; ++resultIt, ++choiceIt, ++rowGroupIt) { - ValueType currentValue = summand ? *summandIt : storm::utility::zero(); - - for (auto elementIte = this->begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { - currentValue += elementIt->getValue() * x[elementIt->getColumn()]; - } - if (choices) { - *choicesIt = 0; - } - - ++rowIt; - ++summandIt; - - for (; *rowIt < *(rowGroupIt + 1); ++rowIt) { - ValueType newValue = summand ? *summandIt : storm::utility::zero(); - for (auto elementIte = this->begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { - newValue += elementIt->getValue() * x[elementIt->getColumn()]; - } - - if ((dir == OptimizationDirection::Minimize && newValue < currentValue) || (dir == OptimizationDirection::Maximize && newValue > currentValue)) { - currentValue = newValue; - if (choices) { - *choiceIt = std::distance(rowIndications.begin(), rowIt) - *rowGroupIt; - } - } - } - - // Finally write value to target vector. - *resultIt = currentValue; - } - }); + tbb::parallel_for(tbb::blocked_range(0, rowGroupIndices.size() - 1, 10), TbbMultAddReduceFunctor(dir, rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices)); } + +#ifdef STORM_HAVE_CARL + template<> + void SparseMatrix::multiplyAndReduceParallel(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); + } +#endif #endif template - void SparseMatrix::multiplyAndReduce(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices, bool allowAliasing, MultiplicationDirection const& multiplicationDirection) const { + void SparseMatrix::multiplyAndReduce(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { - // If the vector and the result are aliases and this is not set to be allowed, we need and temporary vector. + // If the vector and the result are aliases, we need and temporary vector. std::vector* target; std::vector temporary; - bool vectorsAliased = &vector == &result; - if (!allowAliasing && vectorsAliased) { + if (&vector == &result) { STORM_LOG_WARN("Vectors are aliased but are not allowed to be. Using temporary, which is potentially slow."); temporary = std::vector(vector.size()); target = &temporary; - STORM_LOG_WARN_COND(multiplicationDirection != MultiplicationDirection::DontCare, "Not specifying multiplication direction for aliased vectors may yield unexpected results."); } else { target = &result; } - - STORM_LOG_WARN_COND(vectorsAliased || multiplicationDirection == MultiplicationDirection::DontCare, "Setting a multiplication direction for unaliased vectors. Check whether this is intended."); -#ifdef STORM_HAVE_INTELTBB - bool useParallel = !vectorsAliased && multiplicationDirection == MultiplicationDirection::DontCare && this->getNonzeroEntryCount() > 10000; - if (useParallel) { - multiplyAndReduceParallel(dir, rowGroupIndices, vector, summand, *target, choices); - } else { -#endif - if (multiplicationDirection == MultiplicationDirection::Forward || (multiplicationDirection == MultiplicationDirection::DontCare && !vectorsAliased)) { - multiplyAndReduceForward(dir, rowGroupIndices, vector, summand, *target, choices); - } else { - multiplyAndReduceBackward(dir, rowGroupIndices, vector, summand, *target, choices); - } -#ifdef STORM_HAVE_INTELTBB - } -#endif + this->multiplyAndReduceForward(dir, rowGroupIndices, vector, summand, *target, choices); + if (target == &temporary) { std::swap(temporary, result); } diff --git a/src/storm/storage/SparseMatrix.h b/src/storm/storage/SparseMatrix.h index bea7602d8..6cc164567 100644 --- a/src/storm/storage/SparseMatrix.h +++ b/src/storm/storage/SparseMatrix.h @@ -776,23 +776,21 @@ namespace storm { template std::vector getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; - enum class MultiplicationDirection { - Forward, Backward, DontCare - }; - /*! * Multiplies the matrix with the given vector and writes the result to the given result vector. * * @param vector The vector with which to multiply the matrix. * @param result The vector that is supposed to hold the result of the multiplication after the operation. * @param summand If given, this summand will be added to the result of the multiplication. - * @param allowAliasing If set, the vector and result vector may be identical in which case the multiplication - * reuses the updated information in the multiplication (like gauss-seidel). - * @param multiplicationDirection The direction in which to perform the multiplication. If the vector and the - * result vector are aliased, the direction will make a difference as other values will be reused. * @return The product of the matrix and the given vector as the content of the given result vector. */ - void multiplyWithVector(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr, bool allowAliasing = false, MultiplicationDirection const& multiplicationDirection = MultiplicationDirection::DontCare) const; + void multiplyWithVector(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; + + void multiplyWithVectorForward(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; + void multiplyWithVectorBackward(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; +#ifdef STORM_HAVE_INTELTBB + void multiplyWithVectorParallel(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; +#endif /*! * Multiplies the matrix with the given vector, reduces it according to the given direction and and writes @@ -804,14 +802,16 @@ namespace storm { * @param summand If given, this summand will be added to the result of the multiplication. * @param result The vector that is supposed to hold the result of the multiplication after the operation. * @param choices If given, the choices made in the reduction process will be written to this vector. - * @param allowAliasing If set, the vector and result vector may be identical in which case the multiplication - * reuses the updated information in the multiplication (like gauss-seidel). - * @param multiplicationDirection The direction in which to perform the multiplication. If the vector and the - * result vector are aliased, the direction will make a difference as other values will be reused. - * @return The product of the matrix and the given vector as the content of the given result vector. + * @return The resulting vector the content of the given result vector. */ - void multiplyAndReduce(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices, bool allowAliasing = false, MultiplicationDirection const& multiplicationDirection = MultiplicationDirection::DontCare) const; + void multiplyAndReduce(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const; + void multiplyAndReduceForward(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const; + void multiplyAndReduceBackward(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const; +#ifdef STORM_HAVE_INTELTBB + void multiplyAndReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const; +#endif + /*! * Multiplies a single row of the matrix with the given vector and returns the result * @@ -1074,18 +1074,6 @@ namespace storm { */ SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector const& rowGroupIndices, bool insertDiagonalEntries = false) const; - void multiplyWithVectorForward(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; - void multiplyWithVectorBackward(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; -#ifdef STORM_HAVE_INTELTBB - void multiplyWithVectorParallel(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; -#endif - - void multiplyAndReduceForward(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const; - void multiplyAndReduceBackward(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const; -#ifdef STORM_HAVE_INTELTBB - void multiplyAndReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* b, std::vector& result, std::vector* choices) const; -#endif - // The number of rows of the matrix. index_type rowCount; diff --git a/src/storm/utility/VectorHelper.cpp b/src/storm/utility/VectorHelper.cpp new file mode 100644 index 000000000..ac154ca04 --- /dev/null +++ b/src/storm/utility/VectorHelper.cpp @@ -0,0 +1,57 @@ +#include "storm/utility/VectorHelper.h" + +#include "storm/settings/SettingsManager.h" +#include "storm/settings/modules/CoreSettings.h" + +#include "storm/adapters/RationalNumberAdapter.h" +#include "storm/adapters/RationalFunctionAdapter.h" + +#include "storm-config.h" + +#include "storm/utility/vector.h" + +#include "storm/utility/macros.h" +#include "storm/exceptions/InvalidSettingsException.h" +#include "storm/exceptions/NotSupportedException.h" + +namespace storm { + namespace utility { + + template + VectorHelper::VectorHelper() : doParallelize(storm::settings::getModule().isUseIntelTbbSet()) { +#ifndef STORM_HAVE_INTELTBB + STORM_LOG_THROW(!parallelize, storm::exceptions::InvalidSettingsException, "Cannot parallelize without TBB."); +#endif + } + + template + bool VectorHelper::parallelize() const { + return doParallelize; + } + + template + void VectorHelper::reduceVector(storm::solver::OptimizationDirection dir, std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices) const { +#ifdef STORM_HAVE_INTELTBB + if (this->parallelize()) { + storm::utility::vector::reduceVectorMinOrMaxParallel(dir, source, target, rowGrouping, choices); + } else { + storm::utility::vector::reduceVectorMinOrMax(dir, source, target, rowGrouping, choices); + } +#else + storm::utility::vector::reduceVectorMinOrMax(dir, source, target, rowGrouping, choices); +#endif + } + + template<> + void VectorHelper::reduceVector(storm::solver::OptimizationDirection dir, std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices) const { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); + } + + template class VectorHelper; + +#ifdef STORM_HAVE_CARL + template class VectorHelper; + template class VectorHelper; +#endif + } +} diff --git a/src/storm/utility/VectorHelper.h b/src/storm/utility/VectorHelper.h new file mode 100644 index 000000000..293fbaef2 --- /dev/null +++ b/src/storm/utility/VectorHelper.h @@ -0,0 +1,25 @@ +#pragma once + +#include +#include + +#include "storm/solver/OptimizationDirection.h" + +namespace storm { + namespace utility { + + template + class VectorHelper { + public: + VectorHelper(); + + void reduceVector(storm::solver::OptimizationDirection dir, std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices = nullptr) const; + + bool parallelize() const; + + private: + // A flag that stores whether the parallelization to be used. + bool doParallelize; + }; + } +} diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index b894ba4bf..88860c072 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -309,7 +309,21 @@ namespace storm { */ template void applyPointwiseTernary(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, Operation f = Operation()) { + auto firstIt = firstOperand.begin(); + auto firstIte = firstOperand.end(); + auto secondIt = secondOperand.begin(); + auto targetIt = target.begin(); + while (firstIt != firstIte) { + *targetIt = f(*firstIt, *secondIt, *targetIt); + ++targetIt; + ++firstIt; + ++secondIt; + } + } + #ifdef STORM_HAVE_INTELTBB + template + void applyPointwiseTernaryParallel(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, Operation f = Operation()) { tbb::parallel_for(tbb::blocked_range(0, target.size()), [&](tbb::blocked_range const& range) { auto firstIt = firstOperand.begin() + range.begin(); @@ -323,19 +337,8 @@ namespace storm { ++secondIt; } }); -#else - auto firstIt = firstOperand.begin(); - auto firstIte = firstOperand.end(); - auto secondIt = secondOperand.begin(); - auto targetIt = target.begin(); - while (firstIt != firstIte) { - *targetIt = f(*firstIt, *secondIt, *targetIt); - ++targetIt; - ++firstIt; - ++secondIt; - } -#endif } +#endif /*! * Applies the given operation pointwise on the two given vectors and writes the result to the third vector. @@ -347,15 +350,19 @@ namespace storm { */ template void applyPointwise(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, Operation f = Operation()) { + std::transform(firstOperand.begin(), firstOperand.end(), secondOperand.begin(), target.begin(), f); + } + #ifdef STORM_HAVE_INTELTBB + template + void applyPointwiseParallel(std::vector const& firstOperand, std::vector const& secondOperand, std::vector& target, Operation f = Operation()) { tbb::parallel_for(tbb::blocked_range(0, target.size()), [&](tbb::blocked_range const& range) { std::transform(firstOperand.begin() + range.begin(), firstOperand.begin() + range.end(), secondOperand.begin() + range.begin(), target.begin() + range.begin(), f); }); -#else - std::transform(firstOperand.begin(), firstOperand.end(), secondOperand.begin(), target.begin(), f); -#endif } +#endif + /*! * Applies the given function pointwise on the given vector. @@ -366,15 +373,18 @@ namespace storm { */ template void applyPointwise(std::vector const& operand, std::vector& target, Operation f = Operation()) { + std::transform(operand.begin(), operand.end(), target.begin(), f); + } + #ifdef STORM_HAVE_INTELTBB + template + void applyPointwiseParallel(std::vector const& operand, std::vector& target, Operation f = Operation()) { tbb::parallel_for(tbb::blocked_range(0, target.size()), [&](tbb::blocked_range const& range) { std::transform(operand.begin() + range.begin(), operand.begin() + range.end(), target.begin() + range.begin(), f); }); -#else - std::transform(operand.begin(), operand.end(), target.begin(), f); -#endif } +#endif /*! * Adds the two given vectors and writes the result to the target vector. @@ -592,6 +602,73 @@ namespace storm { } return current; } + +#ifdef STORM_HAVE_INTELTBB + template + class TbbReduceVectorFunctor { + public: + TbbReduceVectorFunctor(std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices, Filter const& f) : source(source), target(target), rowGrouping(rowGrouping), choices(choices), f(f) { + // Intentionally left empty. + } + + void operator()(tbb::blocked_range const& range) const { + uint_fast64_t startRow = range.begin(); + uint_fast64_t endRow = range.end(); + + typename std::vector::iterator targetIt = target.begin() + startRow; + typename std::vector::iterator targetIte = target.begin() + endRow; + typename std::vector::const_iterator rowGroupingIt = rowGrouping.begin() + startRow; + typename std::vector::const_iterator sourceIt = source.begin() + *rowGroupingIt; + typename std::vector::const_iterator sourceIte; + typename std::vector::iterator choiceIt; + uint_fast64_t localChoice; + if (choices != nullptr) { + choiceIt = choices->begin() + startRow; + } + + for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt) { + // Only do work if the row group is not empty. + if (*rowGroupingIt != *(rowGroupingIt + 1)) { + *targetIt = *sourceIt; + ++sourceIt; + localChoice = 1; + if (choices != nullptr) { + *choiceIt = 0; + } + + for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { + if (f(*sourceIt, *targetIt)) { + *targetIt = *sourceIt; + if (choices != nullptr) { + *choiceIt = localChoice; + } + } + } + + if (choices != nullptr) { + ++choiceIt; + } + } else { + // Compensate for the 'wrong' move forward in the loop header. + --targetIt; + + // Record dummy choice. + if (choices != nullptr) { + *choiceIt = 0; + ++choiceIt; + } + } + } + } + + private: + std::vector const& source; + std::vector& target; + std::vector const& rowGrouping; + std::vector* choices; + Filter const& f; + }; +#endif /*! * Reduces the given source vector by selecting an element according to the given filter out of each row group. @@ -606,58 +683,6 @@ namespace storm { template void reduceVector(std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices) { Filter f; -#ifdef STORM_HAVE_INTELTBB - tbb::parallel_for(tbb::blocked_range(0, target.size()), - [&](tbb::blocked_range const& range) { - uint_fast64_t startRow = range.begin(); - uint_fast64_t endRow = range.end(); - - typename std::vector::iterator targetIt = target.begin() + startRow; - typename std::vector::iterator targetIte = target.begin() + endRow; - typename std::vector::const_iterator rowGroupingIt = rowGrouping.begin() + startRow; - typename std::vector::const_iterator sourceIt = source.begin() + *rowGroupingIt; - typename std::vector::const_iterator sourceIte; - typename std::vector::iterator choiceIt; - uint_fast64_t localChoice; - if (choices != nullptr) { - choiceIt = choices->begin(); - } - - for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt) { - // Only do work if the row group is not empty. - if (*rowGroupingIt != *(rowGroupingIt + 1)) { - *targetIt = *sourceIt; - ++sourceIt; - localChoice = 1; - if (choices != nullptr) { - *choiceIt = 0; - } - - for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { - if (f(*sourceIt, *targetIt)) { - *targetIt = *sourceIt; - if (choices != nullptr) { - *choiceIt = localChoice; - } - } - } - - if (choices != nullptr) { - ++choiceIt; - } - } else { - // Compensate for the 'wrong' move forward in the loop header. - --targetIt; - - // Record dummy choice. - if (choices != nullptr) { - *choiceIt = 0; - ++choiceIt; - } - } - } - }); -#else typename std::vector::iterator targetIt = target.begin(); typename std::vector::iterator targetIte = target.end(); typename std::vector::const_iterator rowGroupingIt = rowGrouping.begin(); @@ -700,10 +725,14 @@ namespace storm { } } } -#endif } - +#ifdef STORM_HAVE_INTELTBB + template + void reduceVectorParallel(std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices) { + tbb::parallel_for(tbb::blocked_range(0, target.size()), TbbReduceVectorFunctor(source, target, rowGrouping, choices, Filter())); + } +#endif /*! * Reduces the given source vector by selecting the smallest element out of each row group. @@ -718,6 +747,13 @@ namespace storm { reduceVector>(source, target, rowGrouping, choices); } +#ifdef STORM_HAVE_INTELTBB + template + void reduceVectorMinParallel(std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices = nullptr) { + reduceVector>(source, target, rowGrouping, choices); + } +#endif + /*! * Reduces the given source vector by selecting the largest element out of each row group. * @@ -731,6 +767,13 @@ namespace storm { reduceVector>(source, target, rowGrouping, choices); } +#ifdef STORM_HAVE_INTELTBB + template + void reduceVectorMaxParallel(std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices = nullptr) { + reduceVector>(source, target, rowGrouping, choices); + } +#endif + /*! * Reduces the given source vector by selecting either the smallest or the largest out of each row group. * @@ -749,6 +792,17 @@ namespace storm { } } +#ifdef STORM_HAVE_INTELTBB + template + void reduceVectorMinOrMaxParallel(storm::solver::OptimizationDirection dir, std::vector const& source, std::vector& target, std::vector const& rowGrouping, std::vector* choices = nullptr) { + if(dir == storm::solver::OptimizationDirection::Minimize) { + reduceVectorMinParallel(source, target, rowGrouping, choices); + } else { + reduceVectorMaxParallel(source, target, rowGrouping, choices); + } + } +#endif + /*! * Compares the given elements and determines whether they are equal modulo the given precision. The provided flag * additionaly specifies whether the error is computed in relative or absolute terms.