From 2ca83f5f31155b0c359cbabda55c6590aefc25c2 Mon Sep 17 00:00:00 2001 From: dehnert Date: Mon, 3 Dec 2012 15:23:21 +0100 Subject: [PATCH] Added functionality to rapidly extract sub-matrix from our sparse matrix format. --- CMakeLists.txt | 35 ++----- src/mrmc.cpp | 6 +- src/solver/GraphAnalyzer.h | 6 +- src/storage/BitVector.h | 73 ++++++++----- src/storage/SquareSparseMatrix.h | 175 ++++++++++++++++++++++++++++++- 5 files changed, 237 insertions(+), 58 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d83136660..52cf0f1a3 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/src/mrmc.cpp b/src/mrmc.cpp index db2275df6..cabac0a54 100644 --- a/src/mrmc.cpp +++ b/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(new log4cplus::PatternLayout("%-5p - %D{%H:%M:%S} - %b:%L : %m%n"))); + consoleLogAppender->setLayout(std::auto_ptr(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; } diff --git a/src/solver/GraphAnalyzer.h b/src/solver/GraphAnalyzer.h index efe1da1e1..6221250d5 100644 --- a/src/solver/GraphAnalyzer.h +++ b/src/solver/GraphAnalyzer.h @@ -69,7 +69,7 @@ public: * @param universalReachabilityStates The result of the search. */ template - static void getUniversalReachabilityStates(mrmc::models::Dtmc& 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& 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 - static void getUniversalReachabilityStates(mrmc::models::Dtmc& model, const mrmc::storage::BitVector& phiStates, const mrmc::storage::BitVector& psiStates, mrmc::storage::BitVector& alwaysPhiUntilPsiStates) { + static void getAlwaysPhiUntilPsiStates(mrmc::models::Dtmc& 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); } }; diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h index d5fb56c67..79e024c4c 100644 --- a/src/storage/BitVector.h +++ b/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; }; diff --git a/src/storage/SquareSparseMatrix.h b/src/storage/SquareSparseMatrix.h index b83b91eff..ac5260599 100644 --- a/src/storage/SquareSparseMatrix.h +++ b/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(rowCount); Eigen::SparseMatrix* mat = new Eigen::SparseMatrix(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* toGMMXXSparseMatrix() { + // Prepare the resulting matrix. + gmm::csr_matrix* result = new gmm::csr_matrix(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.