Implemented TBB Parallelization Support into SparseMatrix.h

Re-factored Includes in CMake for TBB


Former-commit-id: b5ebf4153a

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

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 <new> // This fixes a potential dependency ordering problem between GMM and TBB
 #	include "tbb/tbb.h"
 #	include <iterator>
@@ -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<typename std::iterator_traits<IT1>::value_type,
 	typename linalg_traits<V>::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<IT1, V> tbbHelper(&v);
 	  tbb::parallel_reduce(forward_range<IT1>(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<L3>::iterator it=vect_begin(l3), ite=vect_end(l3);
     typename linalg_traits<L1>::const_row_iterator
       itr = mat_row_const_begin(l1); 
-#ifdef GMM_USE_TBB
+#ifdef STORM_USE_TBB
     tbb::parallel_for(forward_range_mult<typename linalg_traits<L3>::iterator, typename linalg_traits<L1>::const_row_iterator>(it, ite, itr), tbbHelper_mult_by_row<L1, L2, L3>(&l2));
 #else
 	for (; it != ite; ++it, ++itr)

src/storage/SparseMatrix.h

@@ -9,6 +9,12 @@
 #include <set>
 #include "boost/integer/integer_mask.hpp"
 
+#ifdef STORM_USE_TBB
+#	include <new> // This fixes a potential dependency ordering problem between GMM and TBB
+#	include "tbb/tbb.h"
+#	include <iterator>
+#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<T> const& matrix, uint_fast64_t row = 0) : matrix(matrix), posIndex(matrix.rowIndications[row]), rowIndex(row) {
+		ConstRowsIterator(SparseMatrix<T> 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<T> const& matrix;
+		SparseMatrix<T> 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<T> const& vector, std::vector<T>& result) const {
+#ifdef STORM_USE_TBB
+		tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, result.size()), tbbHelper_MatrixRowVectorScalarProduct<storm::storage::SparseMatrix<T>, std::vector<T>, 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 <typename M, typename V, typename T>
+	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<uint_fast64_t>& 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<T>();
+			
+				// 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
 

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<double> matrix(10, 10);
+	ASSERT_NO_THROW(matrix.initialize(100));
+	double values[100];
+	std::vector<double> 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<double>::MatrixStatus::ReadReady);
+
+	std::vector<double> 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<double> matrix(10, 10);
+	ASSERT_NO_THROW(matrix.initialize(50));
+	double values[100];
+	std::vector<double> 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<double>::MatrixStatus::ReadReady);
+
+	std::vector<double> 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