Added functionality to rapidly extract sub-matrix from our sparse matrix format.

CMakeLists.txt

@@ -25,11 +25,14 @@ message(STATUS "LOG4CPLUS_LIBRARY is ${LOG4CPLUS_LIBRARY}")
 set(EIGEN3_INCLUDE_DIR ${PROJECT_SOURCE_DIR}/resources/3rdparty/eigen)
 message(STATUS "EIGEN3_INCLUDE_DIR is ${EIGEN3_INCLUDE_DIR}")
 
+# Set all Eigen references the version in the repository
+set(GMMXX_INCLUDE_DIR ${PROJECT_SOURCE_DIR}/resources/3rdparty/gmm-4.2/include)
+message(STATUS "GMMXX_INCLUDE_DIR is ${GMMXX_INCLUDE_DIR}")
+
 # Now define all available custom options
 option(DEBUG "Sets whether the DEBUG mode is used" ON)
 option(USE_POPCNT "Sets whether the popcnt instruction is going to be used." ON)
 option(USE_BOOST_STATIC_LIBRARIES "Sets whether the Boost libraries should be linked statically." ON)
-option(USE_BOOST_MT_SUFFIX "Sets whether the the Boost libraries have the suffix mt." OFF)
 
 # If the DEBUG option was turned on, we will target a debug version and a release version otherwise
 if (DEBUG)
@@ -52,7 +55,7 @@ elseif(MSVC)
 
 else(CLANG)
     # Set standard flags for GCC
-    set (CMAKE_CXX_FLAGS "-std=c++11 -stdlib=libc++ -Wall -Werror -pedantic")
+    set (CMAKE_CXX_FLAGS "-std=c++11 -stdlib=libc++ -Wall -Werror -pedantic -Wno-unused-variable")
     
     # Turn on popcnt instruction if desired (yes by default)
     if (USE_POPCNT)
@@ -93,50 +96,34 @@ find_package(Doxygen REQUIRED)
 #set(Boost_USE_STATIC_LIBS        ON)
 set(Boost_USE_MULTITHREADED      ON)
 set(Boost_USE_STATIC_RUNTIME    OFF)
-find_package(Boost REQUIRED)
+find_package(Boost REQUIRED COMPONENTS program_options)
 
 if(Boost_FOUND)
-	if (NOT USE_BOOST_MT_SUFFIX)
-        set(BOOST_PROGRAM_OPTIONS_LIBRARY "boost_program_options")
-	else()
-        set(BOOST_PROGRAM_OPTIONS_LIBRARY "boost_program_options-mt")
-	endif()
-
     if ((NOT Boost_LIBRARY_DIRS) OR ("${Boost_LIBRARY_DIRS}" STREQUAL ""))
         set(Boost_LIBRARY_DIRS "${Boost_INCLUDE_DIRS}/stage/lib")
     endif ()
 	
     message(STATUS "BOOST_INCLUDE_DIRS is ${Boost_INCLUDE_DIRS}")
     message(STATUS "BOOST_LIBRARY_DIRS is ${Boost_LIBRARY_DIRS}")
-	message(STATUS "BOOST_PROGRAM_OPTIONS_LIBRARY is ${BOOST_PROGRAM_OPTIONS_LIBRARY}")
 	
     include_directories(${Boost_INCLUDE_DIRS})
 	link_directories(${Boost_LIBRARY_DIRS})
-	
-	# Prepare the static boost libraries in case they were requested
-    if (USE_BOOST_STATIC_LIBRARIES)
-        add_library(boost_program_options STATIC IMPORTED)
-        set_property(TARGET boost_program_options PROPERTY IMPORTED_LOCATION "${Boost_LIBRARY_DIRS}/lib${BOOST_PROGRAM_OPTIONS_LIBRARY}.a")
-        set_target_properties(boost_program_options PROPERTIES LINKER_LANGUAGE CXX)
-    endif(USE_BOOST_STATIC_LIBRARIES)
 endif(Boost_FOUND)
 
 # Add Eigen to the included directories
 include_directories(${EIGEN3_INCLUDE_DIR})
 
+# Add GMMXX to the included directories
+include_directories(${GMMXX_INCLUDE_DIR})
+
 # Add the executables
 # Must be created *after* Boost was added because of LINK_DIRECTORIES
 add_executable(mrmc ${MRMC_SOURCES} ${MRMC_HEADERS})
 add_executable(mrmc-tests ${MRMC_TEST_SOURCES} ${MRMC_TEST_HEADERS})
 
 # Add target link deps for Boost program options
-if (USE_BOOST_STATIC_LIBRARIES)
-    target_link_libraries(mrmc ${BOOST_PROGRAM_OPTIONS_LIBRARY})
-    target_link_libraries(mrmc-tests ${BOOST_PROGRAM_OPTIONS_LIBRARY})
-else()
-    target_link_libraries(mrmc boost_program_options)
-    target_link_libraries(mrmc-tests boost_program_options)
-endif(USE_BOOST_STATIC_LIBRARIES)
+target_link_libraries(mrmc ${Boost_LIBRARIES})   
+target_link_libraries(mrmc-tests ${Boost_LIBRARIES})
 
 # Add a target to generate API documentation with Doxygen
 if(DOXYGEN_FOUND)

src/mrmc.cpp

@@ -25,6 +25,7 @@
 #include "src/solver/GraphAnalyzer.h"
 #include "src/utility/settings.h"
 #include "Eigen/Sparse"
+#include "gmm/gmm_matrix.h"
 
 #include "log4cplus/logger.h"
 #include "log4cplus/loggingmacros.h"
@@ -49,7 +50,7 @@ void setUpFileLogging() {
 void setUpConsoleLogging() {
 	log4cplus::SharedAppenderPtr consoleLogAppender(new log4cplus::ConsoleAppender());
 	consoleLogAppender->setName("mainConsoleAppender");
-	consoleLogAppender->setLayout(std::auto_ptr<log4cplus::Layout>(new log4cplus::PatternLayout("%-5p - %D{%H:%M:%S} - %b:%L : %m%n")));
+	consoleLogAppender->setLayout(std::auto_ptr<log4cplus::Layout>(new log4cplus::PatternLayout("%-5p - %D{%H:%M:%S} (%r ms) - %b:%L : %m%n")));
 	logger.addAppender(consoleLogAppender);
 }
 
@@ -97,9 +98,6 @@ int main(const int argc, const char* argv[]) {
 
 	dtmc.printModelInformationToStream(std::cout);
 
-	// mrmc::storage::BitVector AU(dtmc.getNumberOfStates());
-	// mrmc::solver::GraphAnalyzer::getUniversalReachabilityStates(dtmc, mrmc::storage::BitVector(dtmc.getNumberOfStates(), true), *dtmc.getLabeledStates("elected"), AU);
-
 	if (s != nullptr) {
 		delete s;
 	}

src/solver/GraphAnalyzer.h

@@ -69,7 +69,7 @@ public:
 	 * @param universalReachabilityStates The result of the search.
 	 */
 	template <class T>
-	static void getUniversalReachabilityStates(mrmc::models::Dtmc<T>& model, const mrmc::storage::BitVector& phiStates, const mrmc::storage::BitVector& psiStates, const mrmc::storage::BitVector& existsPhiUntilPsiStates, mrmc::storage::BitVector& alwaysPhiUntilPsiStates) {
+	static void getAlwaysPhiUntilPsiStates(mrmc::models::Dtmc<T>& model, const mrmc::storage::BitVector& phiStates, const mrmc::storage::BitVector& psiStates, const mrmc::storage::BitVector& existsPhiUntilPsiStates, mrmc::storage::BitVector& alwaysPhiUntilPsiStates) {
 		GraphAnalyzer::getExistsPhiUntilPsiStates(model, ~psiStates, ~existsPhiUntilPsiStates, alwaysPhiUntilPsiStates);
 		alwaysPhiUntilPsiStates.complement();
 	}
@@ -84,10 +84,10 @@ public:
 	 * @param universalReachabilityStates The result of the search.
 	 */
 	template <class T>
-	static void getUniversalReachabilityStates(mrmc::models::Dtmc<T>& model, const mrmc::storage::BitVector& phiStates, const mrmc::storage::BitVector& psiStates, mrmc::storage::BitVector& alwaysPhiUntilPsiStates) {
+	static void getAlwaysPhiUntilPsiStates(mrmc::models::Dtmc<T>& model, const mrmc::storage::BitVector& phiStates, const mrmc::storage::BitVector& psiStates, mrmc::storage::BitVector& alwaysPhiUntilPsiStates) {
 		mrmc::storage::BitVector existsPhiUntilPsiStates(model.getNumberOfStates());
 		GraphAnalyzer::getExistsPhiUntilPsiStates(model, phiStates, psiStates, existsPhiUntilPsiStates);
-		GraphAnalyzer::getUniversalReachabilityStates(model, phiStates, psiStates, existsPhiUntilPsiStates, alwaysPhiUntilPsiStates);
+		GraphAnalyzer::getAlwaysPhiUntilPsiStates(model, phiStates, psiStates, existsPhiUntilPsiStates, alwaysPhiUntilPsiStates);
 	}
 
 };

src/storage/BitVector.h

@@ -33,6 +33,7 @@ public:
 	 * alter the bit vector.
 	 */
 	class constIndexIterator {
+		friend class BitVector;
 	public:
 
 		/*!
@@ -53,14 +54,6 @@ public:
 			}
 		}
 
-		/*!
-		 * Constructs a const iterator over the indices of the set bits in the
-		 * given bit vector, stopping at the given end index.
-		 * @param bitVector The bit vector to iterate over.
-		 * @param endIndex The number of elements to iterate over.
-		 */
-//		constIndexIterator(const BitVector& bitVector, uint_fast64_t endIndex) : constIndexIterator(bitVector, 0, endIndex, true) { }
-
 		/*!
 		 * Increases the position of the iterator to the position of the next bit that
 		 * is set to true.
@@ -97,14 +90,6 @@ public:
 		uint_fast64_t endIndex;
 	};
 
-	//! Constructor
-	/*!
-	 * Constructs a bit vector which can hold the given number of bits and
-	 * initializes all bits to false.
-	 * @param length The number of bits the bit vector should be able to hold.
-	 */
-//	BitVector(uint_fast64_t length) : BitVector(length, false) { }
-
 	//! Constructor
 	/*!
 	 * Constructs a bit vector which can hold the given number of bits and
@@ -112,7 +97,7 @@ public:
 	 * @param length The number of bits the bit vector should be able to hold.
 	 * @param initTrue The initial value of the first |length| bits.
 	 */
-	BitVector(uint_fast64_t length, bool initTrue = false) {
+	BitVector(uint_fast64_t length, bool initTrue = false) : endIterator(*this, length, length, false) {
 		// Check whether the given length is valid.
 		if (length == 0) {
 			LOG4CPLUS_ERROR(logger, "Trying to create bit vector of size 0.");
@@ -147,7 +132,7 @@ public:
 	 * Copy Constructor. Performs a deep copy of the given bit vector.
 	 * @param bv A reference to the bit vector to be copied.
 	 */
-	BitVector(const BitVector &bv) : bucketCount(bv.bucketCount), bitCount(bv.bitCount) {
+	BitVector(const BitVector &bv) : bucketCount(bv.bucketCount), bitCount(bv.bitCount), endIterator(*this, bitCount, bitCount, false) {
 		LOG4CPLUS_WARN(logger, "Invoking copy constructor.");
 		bucketArray = new uint64_t[bucketCount];
 		std::copy(bv.bucketArray, bv.bucketArray + bucketCount, bucketArray);
@@ -178,6 +163,7 @@ public:
 		bitCount = bv.bitCount;
 		bucketArray = new uint64_t[bucketCount];
 		std::copy(bv.bucketArray, bv.bucketArray + bucketCount, bucketArray);
+		updateEndIterator();
 		return *this;
 	}
 
@@ -204,6 +190,8 @@ public:
 			tempArray[i] = 0;
 		}
 
+		updateEndIterator();
+
 		// Dispose of the old bit vector and set the new one.
 		delete[] bucketArray;
 		bucketArray = tempArray;
@@ -395,22 +383,46 @@ public:
 	 * @return The number of bits that are set (to one) in this bit vector.
 	 */
 	uint_fast64_t getNumberOfSetBits() {
-		uint_fast64_t set_bits = 0;
-		for (uint_fast64_t i = 0; i < bucketCount; ++i) {
+		return getNumberOfSetBitsBeforeIndex(bucketCount << 6);
+	}
+
+	uint_fast64_t getNumberOfSetBitsBeforeIndex(uint_fast64_t index) {
+		uint_fast64_t result = 0;
+		// First, count all full buckets.
+		uint_fast64_t bucket = index >> 6;
+		for (uint_fast64_t i = 0; i < bucket; ++i) {
 			// Check if we are using g++ or clang++ and, if so, use the built-in function
 #if (defined (__GNUG__) || defined(__clang__))
-			set_bits += __builtin_popcountll(this->bucketArray[i]);
+			result += __builtin_popcountll(this->bucketArray[i]);
 #else
 			uint_fast32_t cnt;
 			uint_fast64_t bitset = this->bucketArray[i];
 			for (cnt = 0; bitset; cnt++) {
 				bitset &= bitset - 1;
 			}
-			set_bits += cnt;
+			result += cnt;
 #endif
 		}
 
-		return set_bits;
+		// Now check if we have to count part of a bucket.
+		uint64_t tmp = index & mod64mask;
+		if (tmp != 0) {
+			tmp = ((1ll << (tmp & mod64mask)) - 1ll);
+			tmp &= bucketArray[bucket];
+			// Check if we are using g++ or clang++ and, if so, use the built-in function
+#if (defined (__GNUG__) || defined(__clang__))
+			result += __builtin_popcountll(tmp);
+#else
+			uint_fast32_t cnt;
+			uint64_t bitset = tmp;
+			for (cnt = 0; bitset; cnt++) {
+				bitset &= bitset - 1;
+			}
+			result += cnt;
+#endif
+		}
+
+		return result;
 	}
 
 	/*!
@@ -438,8 +450,8 @@ public:
 	/*!
 	 * Returns an iterator pointing at the element past the bit vector.
 	 */
-	constIndexIterator end() const {
-		return constIndexIterator(*this, bitCount, bitCount, false);
+	const constIndexIterator& end() const {
+		return endIterator;
 	}
 
 private:
@@ -497,6 +509,14 @@ private:
 		}
 	}
 
+	/*!
+	 * Updates the end iterator to the correct past-the-end position. Needs
+	 * to be called whenever the size of the bit vector changed.
+	 */
+	void updateEndIterator() {
+		endIterator.currentIndex = bitCount;
+	}
+
 	/*! The number of 64-bit buckets we use as internal storage. */
 	uint_fast64_t bucketCount;
 
@@ -506,6 +526,9 @@ private:
 	/*! Array of 64-bit buckets to store the bits. */
 	uint64_t* bucketArray;
 
+	/*! An iterator marking the end of the bit vector. */
+	constIndexIterator endIterator;
+
 	/*! A bit mask that can be used to reduce a modulo operation to a logical "and".  */
 	static const uint_fast64_t mod64mask = (1 << 6) - 1;
 };

src/storage/SquareSparseMatrix.h

@@ -10,9 +10,11 @@
 #include "src/exceptions/invalid_argument.h"
 #include "src/exceptions/out_of_range.h"
 #include "src/exceptions/file_IO_exception.h"
+#include "src/storage/BitVector.h"
 
 #include "src/misc/const_templates.h"
 #include "Eigen/Sparse"
+#include "gmm/gmm_matrix.h"
 
 #include "log4cplus/logger.h"
 #include "log4cplus/loggingmacros.h"
@@ -471,7 +473,7 @@ public:
 			LOG4CPLUS_ERROR(logger, "Trying to convert a matrix that is not in a readable state to an Eigen matrix.");
 			throw mrmc::exceptions::invalid_state("Trying to convert a matrix that is not in a readable state to an Eigen matrix.");
 		} else {
-			// Create a
+			// Create the resulting matrix.
 			int_fast32_t eigenRows = static_cast<int_fast32_t>(rowCount);
 			Eigen::SparseMatrix<T, Eigen::RowMajor, int_fast32_t>* mat = new Eigen::SparseMatrix<T, Eigen::RowMajor, int_fast32_t>(eigenRows, eigenRows);
 
@@ -488,6 +490,7 @@ public:
 			// than an average row, the other solution might be faster.
 			// The desired conversion method may be set by an appropriate define.
 
+#define MRMC_USE_TRIPLETCONVERT
 #			ifdef MRMC_USE_TRIPLETCONVERT
 
 			// FIXME: Wouldn't it be more efficient to add the elements in
@@ -503,10 +506,12 @@ public:
 			// and add the corresponding triplet.
 			uint_fast64_t rowStart;
 			uint_fast64_t rowEnd;
+			uint_fast64_t zeroCount = 0;
 			for (uint_fast64_t row = 0; row <= rowCount; ++row) {
 				rowStart = rowIndications[row];
 				rowEnd = rowIndications[row + 1];
 				while (rowStart < rowEnd) {
+					if (valueStorage[rowStart] == 0) zeroCount++;
 					tripletList.push_back(IntTriplet(row, columnIndications[rowStart], valueStorage[rowStart]));
 					++rowStart;
 				}
@@ -514,7 +519,8 @@ public:
 
 			// Then add the elements on the diagonal.
 			for (uint_fast64_t i = 0; i <= rowCount; ++i) {
-				tripletList.push_back(IntTriplet(i, i, diagonalStorage[i]));
+				if (diagonalStorage[i] == 0) zeroCount++;
+				// tripletList.push_back(IntTriplet(i, i, diagonalStorage[i]));
 			}
 
 			// Let Eigen create a matrix from the given list of triplets.
@@ -530,16 +536,19 @@ public:
 			// them to the matrix individually.
 			uint_fast64_t rowStart;
 			uint_fast64_t rowEnd;
+			uint_fast64_t count = 0;
 			for (uint_fast64_t row = 0; row < rowCount; ++row) {
 				rowStart = rowIndications[row];
 				rowEnd = rowIndications[row + 1];
 
 				// Insert the element on the diagonal.
 				mat->insert(row, row) = diagonalStorage[row];
+				count++;
 
 				// Insert the elements that are not on the diagonal
 				while (rowStart < rowEnd) {
 					mat->insert(row, columnIndications[rowStart]) = valueStorage[rowStart];
+					count++;
 					++rowStart;
 				}
 			}
@@ -556,6 +565,57 @@ public:
 		return nullptr;
 	}
 
+	/*!
+	 * Converts the matrix into a sparse matrix in the GMMXX format.
+	 * @return A pointer to a column-major sparse matrix in GMMXX format.
+	 */
+	gmm::csr_matrix<T>* toGMMXXSparseMatrix() {
+		// Prepare the resulting matrix.
+		gmm::csr_matrix<T>* result = new gmm::csr_matrix<T>(rowCount, rowCount);
+
+		LOG4CPLUS_INFO(logger, "Starting copy1.");
+		// Copy over the row indications as GMMXX uses the very same internal format.
+		result->jc.reserve(rowCount + 1);
+		std::copy(rowIndications, rowIndications + (rowCount + 1), result->jc.begin());
+		LOG4CPLUS_INFO(logger, "Done copy1.");
+
+		// For the column indications and the actual values, we have to gather
+		// the values in a temporary array first, as we have to integrate
+		// the values from the diagonal.
+		uint_fast64_t realNonZeros = getNonZeroEntryCount() + getDiagonalNonZeroEntryCount();
+		uint_fast64_t* tmpColumnIndicationsArray = new uint_fast64_t[realNonZeros];
+		uint_fast64_t* tmpValueArray = new uint_fast64_t[realNonZeros];
+		T zero(0);
+		uint_fast64_t currentPosition = 0;
+		for (uint_fast64_t i = 0; i < rowCount; ++i) {
+			bool includedDiagonal = false;
+			for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) {
+				if (diagonalStorage[i] != zero && !includedDiagonal && columnIndications[j] > i) {
+					includedDiagonal = true;
+					tmpColumnIndicationsArray[currentPosition] = i;
+					tmpValueArray[currentPosition] = diagonalStorage[i];
+					++currentPosition;
+				}
+				tmpColumnIndicationsArray[currentPosition] = columnIndications[j];
+				tmpValueArray[currentPosition] = valueStorage[j];
+			}
+		}
+
+		LOG4CPLUS_INFO(logger, "Starting copy2.");
+		// Now, we can copy the temporary array to the GMMXX format.
+		result->ir.reserve(realNonZeros);
+		std::copy(tmpColumnIndicationsArray, tmpColumnIndicationsArray + realNonZeros, result->ir.begin());
+		delete[] tmpColumnIndicationsArray;
+
+		// And do the same thing with the actual values.
+		result->pr.resize(realNonZeros);
+		std::copy(tmpValueArray, tmpValueArray + realNonZeros, result->pr.begin());
+		delete[] tmpValueArray;
+		LOG4CPLUS_INFO(logger, "Done copy2.");
+
+		return result;
+	}
+
 	/*!
 	 * Returns the number of non-zero entries that are not on the diagonal.
 	 * @returns The number of non-zero entries that are not on the diagonal.
@@ -564,6 +624,19 @@ public:
 		return nonZeroEntryCount;
 	}
 
+	/*!
+	 * Returns the number of non-zero entries on the diagonal.
+	 * @return The number of non-zero entries on the diagonal.
+	 */
+	uint_fast64_t getDiagonalNonZeroEntryCount() const {
+		uint_fast64_t result = 0;
+		T zero(0);
+		for (uint_fast64_t i = 0; i < rowCount; ++i) {
+			if (diagonalStorage[i] != zero) ++result;
+		}
+		return result;
+	}
+
 	/*!
 	 * This function makes the given row absorbing. This means that all
 	 * entries in will be set to 0 and the value 1 will be written
@@ -594,6 +667,104 @@ public:
 		return true;
 	}
 
+	/*
+	 * Computes the sum of the elements in the given row whose column bits
+	 * are set to one on the given constraint.
+	 * @param row The row whose elements to add.
+	 * @param constraint A bit vector that indicates which columns to add.
+	 * @return The sum of the elements in the given row whose column bits
+	 * are set to one on the given constraint.
+	 */
+	T getConstrainedRowSum(const uint_fast64_t row, const mrmc::storage::BitVector& constraint) {
+		T result(0);
+		for (uint_fast64_t i = rowIndications[row]; i < rowIndications[row + 1]; ++i) {
+			if (constraint.get(columnIndications[i])) {
+				result += valueStorage[i];
+			}
+		}
+		return result;
+	}
+
+	/*!
+	 * Computes a vector in which each element is the sum of those elements in the
+	 * corresponding row whose column bits are set to one in the given constraint.
+	 * @param constraint A bit vector that indicates which columns to add.
+	 * @param resultVector A pointer to the resulting vector that has at least
+	 * as many elements as there are bits set to true in the constraint.
+	 */
+	void getConstrainedRowCountVector(const mrmc::storage::BitVector& constraint, T* resultVector) {
+		for (uint_fast64_t row = 0; row < rowCount; ++row) {
+			resultVector[row] = getConstrainedRowSum(row, constraint);
+		}
+	}
+
+	/*!
+	 * Creates a sub-matrix of the current matrix by dropping all rows and
+	 * columns whose bits are not set to one in the given bit vector.
+	 * @param constraint A bit vector indicating which rows and columns to drop.
+	 * @return A pointer to a sparse matrix that is a sub-matrix of the current one.
+	 */
+	SquareSparseMatrix* getSubmatrix(mrmc::storage::BitVector& constraint) {
+		LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix with " << constraint.getNumberOfSetBits() << " rows.");
+
+		// Check for valid constraint.
+		if (constraint.getNumberOfSetBits() == 0) {
+			LOG4CPLUS_ERROR(logger, "Trying to create a sub-matrix of size 0.");
+			throw mrmc::exceptions::invalid_argument("Trying to create a sub-matrix of size 0.");
+		}
+
+		// First, we need to determine the number of non-zero entries of the
+		// sub-matrix.
+		uint_fast64_t subNonZeroEntries = 0;
+		for (auto rowIndex : constraint) {
+			for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) {
+				if (constraint.get(columnIndications[i])) {
+					++subNonZeroEntries;
+				}
+			}
+		}
+
+		LOG4CPLUS_DEBUG(logger, "Done counting non-zeros.");
+
+		// Create and initialize resulting matrix.
+		SquareSparseMatrix* result = new SquareSparseMatrix(constraint.getNumberOfSetBits());
+		result->initialize(subNonZeroEntries);
+
+		// Create a temporary array that stores for each index whose bit is set
+		// to true the number of bits that were set before that particular index.
+		uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[rowCount];
+		uint_fast64_t lastIndex = 0;
+		uint_fast64_t currentNumberOfSetBits = 0;
+		for (auto index : constraint) {
+			while (lastIndex <= index) {
+				bitsSetBeforeIndex[lastIndex++] = currentNumberOfSetBits;
+			}
+			++currentNumberOfSetBits;
+		}
+
+		// Copy over selected entries.
+		uint_fast64_t rowCount = 0;
+		for (auto rowIndex : constraint) {
+			result->addNextValue(rowCount, rowCount, diagonalStorage[rowIndex]);
+
+			for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) {
+				if (constraint.get(columnIndications[i])) {
+					result->addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[i]], valueStorage[i]);
+				}
+			}
+
+			++rowCount;
+		}
+
+		// Dispose of the temporary array.
+		delete[] bitsSetBeforeIndex;
+
+		// Finalize sub-matrix and return result.
+		result->finalize();
+		LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix.");
+		return result;
+	}
+
 	/*!
 	 * Returns the size of the matrix in memory measured in bytes.
 	 * @return The size of the matrix in memory measured in bytes.