diff --git a/CMakeLists.txt b/CMakeLists.txt index 874f2133c..8afd0b1b0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -215,11 +215,11 @@ if (TBB_FOUND) message(STATUS "Found TBB with Interface Version ${TBB_INTERFACE_VERSION}") if(USE_INTELTBB) - add_definitions(-DGMM_USE_TBB) + add_definitions(-DSTORM_USE_TBB) + + include_directories(${TBB_INCLUDE_DIRS}) + link_directories(${TBB_LIBRARY_DIRS}) endif() - - include_directories(${TBB_INCLUDE_DIRS}) - link_directories(${TBB_LIBRARY_DIRS}) endif(TBB_FOUND) # Add Eigen to the included directories diff --git a/resources/3rdparty/gmm-4.2/include/gmm/gmm_blas.h b/resources/3rdparty/gmm-4.2/include/gmm/gmm_blas.h index 249ae5efb..841a06e75 100644 --- a/resources/3rdparty/gmm-4.2/include/gmm/gmm_blas.h +++ b/resources/3rdparty/gmm-4.2/include/gmm/gmm_blas.h @@ -38,7 +38,7 @@ #ifndef GMM_BLAS_H__ #define GMM_BLAS_H__ -#ifdef GMM_USE_TBB +#ifdef STORM_USE_TBB # include // This fixes a potential dependency ordering problem between GMM and TBB # include "tbb/tbb.h" # include @@ -401,7 +401,7 @@ namespace gmm { return res; } -#ifdef GMM_USE_TBB +#ifdef STORM_USE_TBB /* Official Intel Hint on blocked_range vs. linear iterators: http://software.intel.com/en-us/forums/topic/289505 */ @@ -465,7 +465,7 @@ public: vect_sp_sparse_(IT1 it, IT1 ite, const V &v) { typename strongest_numeric_type::value_type, typename linalg_traits::value_type>::T res(0); -#if defined(GMM_USE_TBB) && defined(GMM_USE_TBB_FOR_INNER) +#if defined(STORM_USE_TBB) && defined(STORM_USE_TBB_FOR_INNER) tbbHelper_vect_sp_sparse tbbHelper(&v); tbb::parallel_reduce(forward_range(it, ite), tbbHelper); res = tbbHelper.my_sum; @@ -1749,7 +1749,7 @@ public: } } -#ifdef GMM_USE_TBB +#ifdef STORM_USE_TBB /* Official Intel Hint on blocked_range vs. linear iterators: http://software.intel.com/en-us/forums/topic/289505 */ @@ -1820,7 +1820,7 @@ public: typename linalg_traits::iterator it=vect_begin(l3), ite=vect_end(l3); typename linalg_traits::const_row_iterator itr = mat_row_const_begin(l1); -#ifdef GMM_USE_TBB +#ifdef STORM_USE_TBB tbb::parallel_for(forward_range_mult::iterator, typename linalg_traits::const_row_iterator>(it, ite, itr), tbbHelper_mult_by_row(&l2)); #else for (; it != ite; ++it, ++itr) diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 8c880806c..f35f172f3 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -9,6 +9,12 @@ #include #include "boost/integer/integer_mask.hpp" +#ifdef STORM_USE_TBB +# include // This fixes a potential dependency ordering problem between GMM and TBB +# include "tbb/tbb.h" +# include +#endif + #include "src/exceptions/InvalidStateException.h" #include "src/exceptions/InvalidArgumentException.h" #include "src/exceptions/OutOfRangeException.h" @@ -76,7 +82,7 @@ public: * * @param matrix The matrix on which this iterator operates. */ - ConstRowsIterator(SparseMatrix const& matrix, uint_fast64_t row = 0) : matrix(matrix), posIndex(matrix.rowIndications[row]), rowIndex(row) { + ConstRowsIterator(SparseMatrix const& matrix, uint_fast64_t row = 0) : matrix(&matrix), posIndex(matrix.rowIndications[row]), rowIndex(row) { // Intentionally left empty. } @@ -87,7 +93,7 @@ public: */ ConstRowsIterator& operator++() { ++posIndex; - if (posIndex >= matrix.rowIndications[rowIndex + 1]) { + if (posIndex >= matrix->rowIndications[rowIndex + 1]) { ++rowIndex; } return *this; @@ -125,6 +131,17 @@ public: return this->posIndex != other.posIndex; } + /*! + * Assignment operator + */ + ConstRowsIterator& operator=( const ConstRowsIterator& rhs) { + this->matrix = rhs.matrix; + this->posIndex = rhs.posIndex; + this->rowIndex = rhs.rowIndex; + + return *this; + } + /*! * Retrieves the row that is associated with the current non-zero element this iterator * points to. @@ -142,7 +159,7 @@ public: * @returns The column of the current non-zero element this iterator points to. */ uint_fast64_t column() const { - return matrix.columnIndications[posIndex]; + return matrix->columnIndications[posIndex]; } /*! @@ -159,7 +176,7 @@ public: * @returns The value of the current non-zero element this iterator points to. */ T const& value() const { - return matrix.valueStorage[posIndex]; + return matrix->valueStorage[posIndex]; } /*! @@ -169,7 +186,7 @@ public: */ void moveToRow(uint_fast64_t row) { this->rowIndex = row; - this->posIndex = matrix.rowIndications[row]; + this->posIndex = matrix->rowIndications[row]; } /*! @@ -179,9 +196,22 @@ public: moveToRow(rowIndex + 1); } + /*! + * Calculates the size of the current row + */ + uint_fast64_t rowSize() { + return (matrix->rowIndications[this->rowIndex + 1] - matrix->rowIndications[this->rowIndex]); + } + + /*! + * Moves the column-pointer forward + */ + void advance(uint_fast64_t count) { + this->posIndex += count; + } private: // A constant reference to the matrix this iterator is associated with. - SparseMatrix const& matrix; + SparseMatrix const* matrix; // The current index in the list of all non-zero elements of the matrix this iterator points to. uint_fast64_t posIndex; @@ -1012,6 +1042,9 @@ public: * vector. */ void multiplyWithVector(std::vector const& vector, std::vector& result) const { +#ifdef STORM_USE_TBB + tbb::parallel_for(tbb::blocked_range(0, result.size()), tbbHelper_MatrixRowVectorScalarProduct, std::vector, T>(this, &vector, &result)); +#else // Initialize two iterators that ConstRowsIterator matrixElementIt(*this); ConstRowsIterator matrixElementIte(*this); @@ -1032,6 +1065,7 @@ public: element += matrixElementIt.value() * vector[matrixElementIt.column()]; } } +#endif } /*! @@ -1375,6 +1409,47 @@ private: } }; +#ifdef STORM_USE_TBB + /*! + * This function is a helper for Parallel Execution of the multipliyWithVector functionality. + * It uses Intels TBB parallel_for paradigm to split up the row/vector multiplication and summation + */ + template + class tbbHelper_MatrixRowVectorScalarProduct { + private: + V * resultVector; + V const* vectorX; + M const* matrixA; + + public: + tbbHelper_MatrixRowVectorScalarProduct(M const* matrixA, V const* vectorX, V * resultVector) : matrixA(matrixA), vectorX(vectorX), resultVector(resultVector) {} + + void operator() (const tbb::blocked_range& r) const { + // Initialize two iterators that + M::ConstRowsIterator matrixElementIt(*matrixA, r.begin()); + M::ConstRowsIterator matrixElementIte(*matrixA, r.begin()); + + for (uint_fast64_t rowNumber = r.begin(); rowNumber != r.end(); ++rowNumber) { + // Put the past-the-end iterator to the correct position. + matrixElementIte.moveToNextRow(); + + // Initialize the result to be 0. + T element = storm::utility::constGetZero(); + + // Perform the scalar product. + for (; matrixElementIt != matrixElementIte; ++matrixElementIt) { + element += matrixElementIt.value() * vectorX->at(matrixElementIt.column()); + } + + // Write back to the result Vector + resultVector->at(rowNumber) = element; + } + } + + }; +#endif + + } // namespace storage } // namespace storm diff --git a/test/functional/storage/SparseMatrixTest.cpp b/test/functional/storage/SparseMatrixTest.cpp index 4ce2a7520..17ab1fffb 100644 --- a/test/functional/storage/SparseMatrixTest.cpp +++ b/test/functional/storage/SparseMatrixTest.cpp @@ -315,3 +315,72 @@ TEST(SparseMatrixTest, ConversionToSparseEigen_RowMajor_SparseMatrixTest) { delete esm; delete ssm; } + +#ifdef STORM_USE_TBB +TEST(SparseMatrixTest, TBBMatrixMultTest) { + storm::storage::SparseMatrix matrix(10, 10); + ASSERT_NO_THROW(matrix.initialize(100)); + double values[100]; + std::vector x; + x.resize(10); + for (uint_fast64_t i = 0; i < 100; ++i) { + values[i] = i + 1.0; + if (i < 10) { + x[i] = 1.0; + } + matrix.addNextValue(i / 10, i % 10, values[i]); + } + ASSERT_NO_THROW(matrix.finalize()); + ASSERT_EQ(matrix.getState(), storm::storage::SparseMatrix::MatrixStatus::ReadReady); + + std::vector result; + result.resize(10); + + matrix.multiplyWithVector(x, result); + + for (uint_fast64_t i = 0; i < 10; ++i) { + double rowSum = 0.0; + for (uint_fast64_t j = 0; j < 10; ++j) { + rowSum += values[10 * i + j]; + } + ASSERT_EQ(rowSum, result.at(i)); + } +} + +TEST(SparseMatrixTest, TBBSparseMatrixMultTest) { + storm::storage::SparseMatrix matrix(10, 10); + ASSERT_NO_THROW(matrix.initialize(50)); + double values[100]; + std::vector x; + x.resize(10); + for (uint_fast64_t i = 0; i < 100; ++i) { + if (i % 2 == 0) { + values[i] = i + 1.0; + matrix.addNextValue(i / 10, i % 10, i + 1.0); + } else { + values[i] = 0.0; + } + + if (i < 10) { + x[i] = 1.0; + } + + } + ASSERT_NO_THROW(matrix.finalize()); + ASSERT_EQ(matrix.getState(), storm::storage::SparseMatrix::MatrixStatus::ReadReady); + + std::vector result; + result.resize(10); + + matrix.multiplyWithVector(x, result); + + for (uint_fast64_t i = 0; i < 10; ++i) { + double rowSum = 0.0; + for (uint_fast64_t j = 0; j < 10; ++j) { + rowSum += values[10 * i + j]; + } + ASSERT_EQ(rowSum, result.at(i)); + } +} + +#endif // STORM_USE_TBB