diff --git a/CMakeLists.txt b/CMakeLists.txt index 06b47443a..8c3ad57ab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,19 +6,24 @@ project (MRMC-cpp CXX C) set (MRMC_CPP_VERSION_MAJOR 1) set (MRMC_CPP_VERSION_MINOR 0) +option(DEBUG "Sets whether the DEBUG mode is used" ON) option(DEFINE_UNIX "Defines the UNIX flag for compilation." OFF) +option(USE_POPCNT "Sets whether the popcnt instruction is going to be used." OFF) + +if (DEBUG) + set (CMAKE_BUILD_TYPE "DEBUG") +else() + set (CMAKE_BUILD_TYPE "RELEASE") +endif() #Configurations for GCC if(CMAKE_COMPILER_IS_GNUCC) - set(MRMC_GCC_NO_DEBUG_SYMBOLS OFF CACHE BOOL "Whether debug symbols should be included in MRMC with GCC, set ON for better performance") - #Using C++11 - set (CMAKE_CXX_FLAGS "-std=c++0x") + set (CMAKE_CXX_FLAGS "-std=c++0x -Wall -pedantic") - #Debug symbols - if (NOT MRMC_GCC_NO_DEBUG_SYMBOLS) - set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g") - message(STATUS "Using Debug Symbols for GCC") - endif() + # Turn on popcnt instruction if possible + if (USE_POPCNT) + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -msse4.2") + endif(USE_POPCNT) endif() # configure a header file to pass some of the CMake settings diff --git a/src/models/atomic_propositions_labeling.h b/src/models/atomic_propositions_labeling.h index fde4c67e2..83b7ced16 100644 --- a/src/models/atomic_propositions_labeling.h +++ b/src/models/atomic_propositions_labeling.h @@ -58,12 +58,8 @@ public: * @param state_count The number of states of the model. * @param ap_count The number of atomic propositions. */ - AtomicPropositionsLabeling(const uint_fast32_t state_count, - const uint_fast32_t ap_count) { - // add one for 1-bound indices - this->state_count = state_count + 1; - this->ap_count = ap_count; - this->aps_current = 0; + AtomicPropositionsLabeling(const uint_fast32_t state_count, const uint_fast32_t ap_count) + : state_count(state_count), ap_count(ap_count), aps_current(0) { this->single_labelings = new SingleAtomicPropositionLabeling*[ap_count]; for (uint_fast32_t i = 0; i < ap_count; ++i) { this->single_labelings[i] = new SingleAtomicPropositionLabeling(state_count); @@ -75,11 +71,11 @@ public: * Copy Constructor. Performs a deep copy of the given atomic proposition * labeling. */ - AtomicPropositionsLabeling(const AtomicPropositionsLabeling& atomic_propositions_labeling) : - state_count(atomic_propositions_labeling.state_count), - ap_count(atomic_propositions_labeling.ap_count), - aps_current(atomic_propositions_labeling.aps_current), - name_to_labeling_map(atomic_propositions_labeling.name_to_labeling_map) { + AtomicPropositionsLabeling(const AtomicPropositionsLabeling& atomic_propositions_labeling) + : state_count(atomic_propositions_labeling.state_count), + ap_count(atomic_propositions_labeling.ap_count), + aps_current(atomic_propositions_labeling.aps_current), + name_to_labeling_map(atomic_propositions_labeling.name_to_labeling_map) { this->single_labelings = new SingleAtomicPropositionLabeling*[ap_count]; for (uint_fast32_t i = 0; i < ap_count; ++i) { this->single_labelings[i] = new SingleAtomicPropositionLabeling( diff --git a/src/models/backward_transitions.h b/src/models/backward_transitions.h index 8a4aca260..a99f0ae3e 100644 --- a/src/models/backward_transitions.h +++ b/src/models/backward_transitions.h @@ -8,7 +8,7 @@ #ifndef BACKWARD_TRANSITIONS_H_ #define BACKWARD_TRANSITIONS_H_ -#include +#include #include "src/sparse/static_sparse_matrix.h" namespace mrmc { @@ -23,63 +23,61 @@ template class BackwardTransitions { public: + /*! + * Just typedef the iterator as a pointer to the index type. + */ + typedef const uint_fast64_t * const state_predecessor_iterator; + //! Constructor /*! * Constructs a backward transitions object from the given sparse matrix * representing the (forward) transition relation. - * @param transition_matrix The (1-based) matrix representing the transition + * @param transition_matrix The (0-based) matrix representing the transition * relation. */ - BackwardTransitions(mrmc::sparse::StaticSparseMatrix* transition_matrix) { - this->state_indices_list = new uint_fast64_t[transition_matrix->getRowCount() + 2]; - this->predecessor_list = new uint_fast64_t[transition_matrix->getNonZeroEntryCount()]; + BackwardTransitions(mrmc::sparse::StaticSparseMatrix* transitionMatrix) + : numberOfStates(transitionMatrix->getRowCount()), + numberOfNonZeroTransitions(transitionMatrix->getNonZeroEntryCount()) { + this->state_indices_list = new uint_fast64_t[numberOfStates + 1]; + this->predecessor_list = new uint_fast64_t[numberOfNonZeroTransitions]; // First, we need to count how many backward transitions each state has. // NOTE: We disregard the diagonal here, as we only consider "true" // predecessors. - // Start by counting all but the last row. - uint_fast64_t* row_indications = transition_matrix->getRowIndicationsPointer(); - uint_fast64_t* column_indications = transition_matrix->getColumnIndicationsPointer(); - for (uint_fast64_t i = 1; i < transition_matrix->getRowCount(); i++) { - for (uint_fast64_t j = row_indications[i]; j < row_indications[i + 1]; j++) { - this->state_indices_list[column_indications[j]]++; + for (uint_fast64_t i = 0; i <= numberOfStates; i++) { + for (auto rowIt = transitionMatrix->beginConstColumnNoDiagIterator(i); rowIt != transitionMatrix->endConstColumnNoDiagIterator(i); ++rowIt) { + this->state_indices_list[*rowIt + 1]++; } } - // Add the last row individually, as the comparison bound in the for-loop - // is different in this case. - for (uint_fast64_t j = row_indications[transition_matrix->getRowCount()]; j < transition_matrix->getNonZeroEntryCount(); j++) { - this->state_indices_list[column_indications[j]]++; - } // Now compute the accumulated offsets. - for (uint_fast64_t i = 1; i < transition_matrix->getRowCount() + 1; i++) { + for (uint_fast64_t i = 1; i <= numberOfStates; i++) { this->state_indices_list[i] = this->state_indices_list[i - 1] + this->state_indices_list[i]; } // Put a sentinel element at the end of the indices list. This way, // for each state i the range of indices can be read off between // state_indices_list[i] and state_indices_list[i + 1]. - this->state_indices_list[transition_matrix->getRowCount() + 1] = this->state_indices_list[transition_matrix->getRowCount()]; + this->state_indices_list[numberOfStates + 1] = numberOfNonZeroTransitions; // Create an array that stores the next index for each state. Initially // this corresponds to the previously computed accumulated offsets. - uint_fast64_t* next_state_index_list = new uint_fast64_t[transition_matrix->getRowCount() + 1]; - memcpy(next_state_index_list, state_indices_list, (transition_matrix->getRowCount() + 1) * sizeof(uint_fast64_t)); + uint_fast64_t* next_state_index_list = new uint_fast64_t[numberOfStates + 1]; + memcpy(next_state_index_list, state_indices_list, (numberOfStates + 1) * sizeof(uint_fast64_t)); // Now we are ready to actually fill in the list of predecessors for // every state. Again, we start by considering all but the last row. - for (uint_fast64_t i = 1; i < transition_matrix->getRowCount(); i++) { - for (uint_fast64_t j = row_indications[i]; j < row_indications[i + 1]; j++) { - this->predecessor_list[next_state_index_list[i]++] = column_indications[j]; + for (uint_fast64_t i = 0; i <= numberOfStates; i++) { + for (auto rowIt = transitionMatrix->beginConstColumnNoDiagIterator(i); rowIt != transitionMatrix->endConstColumnNoDiagIterator(i); ++rowIt) { + this->predecessor_list[next_state_index_list[*rowIt]++] = i; } } - // Add the last row individually, as the comparison bound in the for-loop - // is different in this case. - for (uint_fast64_t j = row_indications[transition_matrix->getRowCount()]; j < transition_matrix->getNonZeroEntryCount(); j++) { - this->state_indices_list[next_state_index_list[transition_matrix->getRowCount()]++] = column_indications[j]; - } } + //! Destructor + /*! + * Destructor. Frees the internal storage. + */ ~BackwardTransitions() { if (this->predecessor_list != NULL) { delete[] this->predecessor_list; @@ -89,26 +87,25 @@ public: } } - class const_iterator : public std::iterator { - public: - const_iterator(uint_fast64_t* ptr) : ptr_(ptr) { } - - const_iterator operator++() { const_iterator i = *this; ptr_++; return i; } - const_iterator operator++(int offset) { ptr_ += offset; return *this; } - const uint_fast64_t& operator*() { return *ptr_; } - const uint_fast64_t* operator->() { return ptr_; } - bool operator==(const const_iterator& rhs) { return ptr_ == rhs.ptr_; } - bool operator!=(const const_iterator& rhs) { return ptr_ != rhs.ptr_; } - private: - uint_fast64_t* ptr_; - }; - - const_iterator beginPredecessorIterator(uint_fast64_t state) const { - return const_iterator(&(this->predecessor_list[this->state_indices_list[state]])); + /*! + * Returns an iterator to the predecessors of the given states. + * @param state The state for which to get the predecessor iterator. + * @return An iterator to the predecessors of the given states. + */ + state_predecessor_iterator beginStatePredecessorIterator(uint_fast64_t state) const { + return this->predecessor_list + this->state_indices_list[state]; } - const_iterator endPredecessorIterator(uint_fast64_t state) const { - return const_iterator(&(this->predecessor_list[this->state_indices_list[state + 1]])); + + /*! + * Returns an iterator referring to the element after the predecessors of + * the given state. + * @param row The state for which to get the iterator. + * @return An iterator referring to the element after the predecessors of + * the given state. + */ + state_predecessor_iterator endStatePredecessorIterator(uint_fast64_t state) const { + return this->predecessor_list + this->state_indices_list[state + 1]; } private: @@ -121,6 +118,17 @@ private: */ uint_fast64_t* state_indices_list; + /*! + * Store the number of states to determine the highest index at which the + * state_indices_list may be accessed. + */ + uint_fast64_t numberOfStates; + + /*! + * Store the number of non-zero transition entries to determine the highest + * index at which the predecessor_list may be accessed. + */ + uint_fast64_t numberOfNonZeroTransitions; }; } // namespace models diff --git a/src/mrmc-cpp.cpp b/src/mrmc-cpp.cpp index d28bcb2d0..e9db4cd67 100644 --- a/src/mrmc-cpp.cpp +++ b/src/mrmc-cpp.cpp @@ -23,7 +23,7 @@ PANTHEIOS_EXTERN_C PAN_CHAR_T const PANTHEIOS_FE_PROCESS_IDENTITY[] = "mrmc-cpp"; #include "MRMCConfig.h" - +#include "src/models/dtmc.h" #include "src/sparse/static_sparse_matrix.h" #include "src/models/atomic_propositions_labeling.h" #include "src/parser/read_lab_file.h" @@ -35,80 +35,16 @@ int main(int argc, char* argv[]) { pantheios_be_file_setFilePath("log.all"); pantheios::log_INFORMATIONAL("MRMC-Cpp started."); - if (argc < 2) { + if (argc < 3) { std::cout << "Required argument #1 inputTraFile.tra not found!" << std::endl; exit(-1); } - mrmc::sparse::StaticSparseMatrix *ssm; - mrmc::sparse::StaticSparseMatrix *ssmBack; - - // 1. Create an instance of the platformstl::performance_counter. (On - // UNIX this will resolve to unixstl::performance_counter; on Win32 it - // will resolve to winstl::performance_counter.) - platformstl::performance_counter counter; - - std::cout << "Building Matrix from File..." << std::endl; - - // 2. Begin the measurement - counter.start(); - - ssm = mrmc::parser::read_tra_file(argv[1]); - //lab = mrmc::parser::read_lab_file(20302, "csl_unbounded_until_sim_06.lab"); - - counter.stop(); - - std::cout << "Loaded " << ssm->getNonZeroEntryCount() << " entries into (" << ssm->getRowCount() << " x " << ssm->getRowCount() << ")." << std::endl; - std::cout << "Time needed was " << (unsigned long long)counter.get_microseconds() << std::endl; - - - - counter.start(); - auto esm = ssm->toEigenSparseMatrix(); - counter.stop(); - - std::cout << "Loaded " << esm->nonZeros() << " entries into (" << esm->rows() << " x " << esm->cols() << ")." << std::endl; - std::cout << "Time needed was " << (unsigned long long)counter.get_microseconds() << std::endl; + mrmc::sparse::StaticSparseMatrix* probMatrix = mrmc::parser::read_tra_file(argv[1]); + mrmc::models::AtomicPropositionsLabeling* labeling = mrmc::parser::read_lab_file(probMatrix->getRowCount(), argv[2]); + mrmc::models::Dtmc dtmc(probMatrix, labeling); - ssmBack = new mrmc::sparse::StaticSparseMatrix(esm->rows()); - - counter.start(); - ssmBack->initialize(*esm); - counter.stop(); - - std::cout << "Loaded " << ssmBack->getNonZeroEntryCount() << " entries into (" << ssmBack->getRowCount() << " x " << ssmBack->getRowCount() << ")." << std::endl; - std::cout << "Time needed was " << (unsigned long long)counter.get_microseconds() << std::endl; - - delete ssm; - delete ssmBack; - - std::cout << std::endl; - std::cout << ":: C-style vs. C++ Style Array init Showdown ::" << std::endl; - - uint_fast64_t *iArray = NULL; - uint_fast32_t i = 0; - uint_fast32_t j = 0; - - counter.start(); - for (i = 0; i < 10000; ++i) { - iArray = (uint_fast64_t*)malloc(2097152 * sizeof(uint_fast64_t)); - for (j = 0; j < 2097152; ++j) { - iArray[j] = 0; - } - free(iArray); - } - counter.stop(); - - std::cout << "C-style: " << (unsigned long long)counter.get_microseconds() << std::endl; - - counter.start(); - for (i = 0; i < 10000; ++i) { - iArray = new uint_fast64_t[2097152](); - delete[] iArray; - } - counter.stop(); - std::cout << "Cpp-style: " << (unsigned long long)counter.get_microseconds() << std::endl; return 0; } diff --git a/src/parser/read_tra_file.cpp b/src/parser/read_tra_file.cpp index 4a62af89c..5d652758b 100644 --- a/src/parser/read_tra_file.cpp +++ b/src/parser/read_tra_file.cpp @@ -135,7 +135,7 @@ sparse::StaticSparseMatrix * read_tra_file(const char * filename) { * Memory for diagonal elements is automatically allocated, hence only the number of non-diagonal * non-zero elements has to be specified (which is non_zero, computed by make_first_pass) */ - sp = new sparse::StaticSparseMatrix(rows); + sp = new sparse::StaticSparseMatrix(rows + 1); if ( NULL == sp ) { throw std::bad_alloc(); return NULL; diff --git a/src/sparse/static_sparse_matrix.h b/src/sparse/static_sparse_matrix.h index e5d586c37..b778fd352 100644 --- a/src/sparse/static_sparse_matrix.h +++ b/src/sparse/static_sparse_matrix.h @@ -21,127 +21,142 @@ namespace mrmc { namespace sparse { -//! A sparse Matrix for DTMCs with a constant number of non-zero entries on the non-diagonal fields. -/*! - The sparse matrix used for calculation on the DTMCs. - Addressing is NOT zero-based! The valid range for getValue and addNextValue is 1 to rows (first argument to constructor). +/*! + * A sparse matrix class with a constant number of non-zero entries on the non-diagonal fields + * and a separate dense storage for the diagonal elements. + * NOTE: Addressing *is* zero-based, so the valid range for getValue and addNextValue is 0..(rows - 1) + * where rows is the first argument to the constructor. */ -template +template class StaticSparseMatrix { - public: - //! Status enum +public: + /*! - An Enum representing the internal state of the Matrix. - After creating the Matrix using the Constructor, the Object is in state UnInitialized. After calling initialize(), that state changes to Initialized and after all entries have been entered and finalize() has been called, to ReadReady. - Should a critical error occur in any of the former functions, the state will change to Error. - @see getState() - @see isReadReady() - @see isInitialized() - @see hasError() + * If we only want to iterate over the columns of the non-zero entries of + * a row, we can simply iterate over the array (part) itself. */ - enum MatrixStatus { - Error = -1, - UnInitialized = 0, - Initialized = 1, - ReadReady = 2, - }; + typedef const uint_fast64_t * const constIndexIterator; + /*! + * An enum representing the internal state of the Matrix. + * After creating the Matrix using the Constructor, the Object is in state UnInitialized. After calling initialize(), that state changes to Initialized and after all entries have been entered and finalize() has been called, to ReadReady. + * Should a critical error occur in any of the former functions, the state will change to Error. + * @see getState() + * @see isReadReady() + * @see isInitialized() + * @see hasError() + */ + enum MatrixStatus { + Error = -1, UnInitialized = 0, Initialized = 1, ReadReady = 2, + }; //! Constructor - /*! - \param rows Row-Count and therefore column-count of the square matrix - */ - StaticSparseMatrix(uint_fast64_t rows) { - // Using direct access instead of setState() because of undefined initialization value - // setState() stays in Error should Error be the current value - internal_status = MatrixStatus::UnInitialized; - current_size = 0; - - value_storage = NULL; - diagonal_storage = NULL; - column_indications = NULL; - row_indications = NULL; - - row_count = rows; - non_zero_entry_count = 0; - - //initialize(rows, non_zero_entries); - } + /*! + * Constructs a sparse matrix object with the given number of rows. + * @param rows The number of rows of the matrix + */ + StaticSparseMatrix(uint_fast64_t rows) + : internalStatus(MatrixStatus::UnInitialized), + currentSize(0), lastRow(0), valueStorage(nullptr), + diagonalStorage(nullptr), columnIndications(nullptr), + rowIndications(nullptr), rowCount(rows), nonZeroEntryCount(0) { } //! Copy Constructor /*! - Copy Constructor. Creates an exact copy of the source sparse matrix ssm. Modification of either matrix does not affect the other. - @param ssm A reference to the matrix that should be copied from + * Copy Constructor. Performs a deep copy of the given sparse matrix. + * @param ssm A reference to the matrix to be copied. */ - StaticSparseMatrix(const StaticSparseMatrix &ssm) : internal_status(ssm.internal_status), current_size(ssm.current_size), row_count(ssm.row_count), non_zero_entry_count(ssm.non_zero_entry_count) - { - pantheios::log_DEBUG("StaticSparseMatrix::CopyCTor: Using Copy() Ctor."); + StaticSparseMatrix(const StaticSparseMatrix &ssm) + : internalStatus(ssm.internalStatus), + currentSize(ssm.currentSize), lastRow(ssm.lastRow), + rowCount(ssm.rowCount), + nonZeroEntryCount(ssm.nonZeroEntryCount) { + pantheios::log_DEBUG("StaticSparseMatrix::CopyConstructor: Using copy constructor."); + // Check whether copying the matrix is safe. if (!ssm.hasError()) { - pantheios::log_ERROR("StaticSparseMatrix::CopyCtor: Throwing invalid_argument: Can not Copy from matrix in Error state."); + pantheios::log_ERROR("StaticSparseMatrix::CopyConstructor: Throwing invalid_argument: Can not copy from matrix in error state."); throw mrmc::exceptions::invalid_argument(); } else { + // Try to prepare the internal storage and throw an error in case + // of a failure. if (!prepareInternalStorage()) { - pantheios::log_ERROR("StaticSparseMatrix::CopyCtor: Throwing bad_alloc: memory allocation failed."); + pantheios::log_ERROR("StaticSparseMatrix::CopyConstructor: Throwing bad_alloc: memory allocation failed."); throw std::bad_alloc(); } else { - for (uint_fast64_t i = 0; i < non_zero_entry_count; ++i) { + // Now that all storages have been prepared, copy over all + // elements. Start by copying the elements of type value and + // copy them seperately in order to invoke copy their copy + // constructor. This may not be necessary, but it is safer to + // do so in any case. + for (uint_fast64_t i = 0; i < nonZeroEntryCount; ++i) { // use T() to force use of the copy constructor for complex T types - value_storage[i] = T(ssm.value_storage[i]); + valueStorage[i] = T(ssm.valueStorage[i]); } - for (uint_fast64_t i = 0; i <= row_count; ++i) { - // same here - diagonal_storage[i] = T(ssm.diagonal_storage[i]); + for (uint_fast64_t i = 0; i <= rowCount; ++i) { + // use T() to force use of the copy constructor for complex T types + diagonalStorage[i] = T(ssm.diagonalStorage[i]); } - memcpy(column_indications, ssm.column_indications, sizeof(column_indications[0]) * non_zero_entry_count); - memcpy(row_indications, ssm.row_indications, sizeof(row_indications[0]) * (row_count + 1)); + // The elements that are not of the value type but rather the + // index type may be copied with memcpy. + memcpy(columnIndications, ssm.columnIndications, sizeof(columnIndications[0]) * nonZeroEntryCount); + memcpy(rowIndications, ssm.rowIndications, sizeof(rowIndications[0]) * (rowCount + 1)); } } } + //! Destructor + /*! + * Destructor. Performs deletion of the reserved storage arrays. + */ ~StaticSparseMatrix() { setState(MatrixStatus::UnInitialized); - if (value_storage != NULL) { + if (valueStorage != NULL) { //free(value_storage); - delete[] value_storage; + delete[] valueStorage; } - if (column_indications != NULL) { + if (columnIndications != NULL) { //free(column_indications); - delete[] column_indications; + delete[] columnIndications; } - if (row_indications != NULL) { + if (rowIndications != NULL) { //free(row_indications); - delete[] row_indications; + delete[] rowIndications; } - if (diagonal_storage != NULL) { + if (diagonalStorage != NULL) { //free(diagonal_storage); - delete[] diagonal_storage; + delete[] diagonalStorage; } } - //! Mandatory initialization of the matrix, variant for initialize(), addNextValue() and finalize() /*! - Mandatory initialization of the matrix, must be called before using any other member function. - This version is to be used together with addNextValue(). - For initialization from a Eigen SparseMatrix, use initialize(Eigen::SparseMatrix &). - \param non_zero_entries The exact count of entries that will be submitted through addNextValue *excluding* those on the diagonal (A_{i,j} with i = j) + * Initializes the sparse matrix with the given number of non-zero entries + * and prepares it for use with addNextValue() and finalize(). + * NOTE: Calling this method before any other member function is mandatory. + * This version is to be used together with addNextValue(). For + * initialization from an Eigen SparseMatrix, use initialize(Eigen::SparseMatrix &). + * @param nonZeroEntries The number of non-zero entries that are not on the + * diagonal. */ - void initialize(uint_fast64_t non_zero_entries) { - if (internal_status != MatrixStatus::UnInitialized) { + void initialize(uint_fast64_t nonZeroEntries) { + // Check whether initializing the matrix is safe. + if (internalStatus != MatrixStatus::UnInitialized) { triggerErrorState(); - pantheios::log_ERROR("StaticSparseMatrix::initialize: Throwing invalid state for status flag != 0 (is ", pantheios::integer(internal_status)," - Already initialized?"); + pantheios::log_ERROR("StaticSparseMatrix::initialize: Throwing invalid state for status flag != 0 (is ", pantheios::integer(internalStatus), " - Already initialized?"); throw mrmc::exceptions::invalid_state("StaticSparseMatrix::initialize: Invalid state for status flag != 0 - Already initialized?"); - } else if (row_count == 0) { + } else if (rowCount == 0) { triggerErrorState(); pantheios::log_ERROR("StaticSparseMatrix::initialize: Throwing invalid_argument for row_count = 0"); throw mrmc::exceptions::invalid_argument("mrmc::StaticSparseMatrix::initialize: Matrix with 0 rows is not reasonable"); - } else if (((row_count * row_count) - row_count) < non_zero_entries) { + } else if (((rowCount * rowCount) - rowCount) < nonZeroEntries) { triggerErrorState(); pantheios::log_ERROR("StaticSparseMatrix::initialize: Throwing invalid_argument: More non-zero entries than entries in target matrix"); throw mrmc::exceptions::invalid_argument("mrmc::StaticSparseMatrix::initialize: More non-zero entries than entries in target matrix"); } else { - non_zero_entry_count = non_zero_entries; - last_row = 0; + // If it is safe, initialize necessary members and prepare the + // internal storage. + nonZeroEntryCount = nonZeroEntries; + lastRow = 0; if (!prepareInternalStorage()) { triggerErrorState(); @@ -153,70 +168,91 @@ class StaticSparseMatrix { } } - //! Mandatory initialization of the matrix, variant for initialize(), addNextValue() and finalize() /*! - Mandatory initialization of the matrix, must be called before using any other member function. - This version is to be used for initialization from a Eigen SparseMatrix, use initialize(uint_fast32_t) for addNextValue. - \param eigen_sparse_matrix The Eigen Sparse Matrix to be copied/ initialized from. MUST BE in compressed form! + * Initializes the sparse matrix with the given Eigen sparse matrix. + * NOTE: Calling this method before any other member function is mandatory. + * This version is only to be used when copying an Eigen sparse matrix. For + * initialization with addNextValue() and finalize() use initialize(uint_fast32_t) + * instead. + * @param eigenSparseMatrix The Eigen sparse matrix to be copied. + * *NOTE* Has to be in compressed form! */ template - void initialize(const Eigen::SparseMatrix &eigen_sparse_matrix) { - if (!eigen_sparse_matrix.isCompressed()) { + void initialize(const Eigen::SparseMatrix& eigenSparseMatrix) { + // Throw an error in case the matrix is not in compressed format. + if (!eigenSparseMatrix.isCompressed()) { triggerErrorState(); pantheios::log_ERROR("StaticSparseMatrix::initialize: Throwing invalid_argument: eigen_sparse_matrix is not in Compressed form."); throw mrmc::exceptions::invalid_argument("StaticSparseMatrix::initialize: Throwing invalid_argument: eigen_sparse_matrix is not in Compressed form."); } - non_zero_entry_count = getEigenSparseMatrixCorrectNonZeroEntryCount(eigen_sparse_matrix); - last_row = 0; + // Compute the actual (i.e. non-diagonal) number of non-zero entries. + nonZeroEntryCount = getEigenSparseMatrixCorrectNonZeroEntryCount(eigenSparseMatrix); + lastRow = 0; + // Try to prepare the internal storage and throw an error in case of + // failure. if (!prepareInternalStorage()) { triggerErrorState(); - pantheios::log_ERROR("StaticSparseMatrix::initialize: Throwing bad_alloc: memory allocation failed"); + pantheios::log_ERROR( + "StaticSparseMatrix::initialize: Throwing bad_alloc: memory allocation failed"); throw std::bad_alloc(); } else { - // easy case, we can simply copy the data - // RowMajor: Easy, ColMajor: Hmm. But how to detect? - const T* valuePtr = eigen_sparse_matrix.valuePtr(); - const _Index* indexPtr = eigen_sparse_matrix.innerIndexPtr(); - const _Index* outerPtr = eigen_sparse_matrix.outerIndexPtr(); - - - const _Index entryCount = eigen_sparse_matrix.nonZeros(); - const _Index outerCount = eigen_sparse_matrix.outerSize(); - if (isEigenRowMajor(eigen_sparse_matrix)) { - // Easy case, all data can be copied with some minor changes. - // We have to separate diagonal entries from others - for (_Index row = 1; row <= outerCount; ++row) { - for (_Index col = outerPtr[row - 1]; col < outerPtr[row]; ++col) { - addNextValue(row, indexPtr[col] + 1, valuePtr[col]); + // Get necessary pointers to the contents of the Eigen matrix. + const T* valuePtr = eigenSparseMatrix.valuePtr(); + const _Index* indexPtr = eigenSparseMatrix.innerIndexPtr(); + const _Index* outerPtr = eigenSparseMatrix.outerIndexPtr(); + + // If the given matrix is in RowMajor format, copying can simply + // be done by adding all values in order. + // Direct copying is, however, prevented because we have to + // separate the diagonal entries from others. + if (isEigenRowMajor(eigenSparseMatrix)) { + // Because of the RowMajor format outerSize evaluates to the + // number of rows. + const _Index rowCount = eigenSparseMatrix.outerSize(); + for (_Index row = 0; row < rowCount; ++row) { + for (_Index col = outerPtr[row]; col < outerPtr[row + 1]; ++col) { + addNextValue(row, indexPtr[col], valuePtr[col]); } } } else { - // temp copies, anyone? - const _Index eigen_col_count = eigen_sparse_matrix.cols(); - const _Index eigen_row_count = eigen_sparse_matrix.rows(); - - // initialise all column-start positions to known lower boundarys - _Index* positions = new _Index[eigen_col_count](); - for (_Index i = 0; i < eigen_col_count; ++i) { + const _Index entryCount = eigenSparseMatrix.nonZeros(); + // Because of the ColMajor format outerSize evaluates to the + // number of columns. + const _Index colCount = eigenSparseMatrix.outerSize(); + + // Create an array to remember which elements have to still + // be searched in each column and initialize it with the starting + // index for every column. + _Index* positions = new _Index[colCount](); + for (_Index i = 0; i < colCount; ++i) { positions[i] = outerPtr[i]; } + // Now copy the elements. As the matrix is in ColMajor format, + // we need to iterate over the columns to find the next non-zero + // entry. int i = 0; int currentRow = 0; int currentColumn = 0; while (i < entryCount) { - if ((positions[currentColumn] < outerPtr[currentColumn + 1]) && (indexPtr[positions[currentColumn]] == currentRow)) { - addNextValue(currentRow + 1, currentColumn + 1, valuePtr[positions[currentColumn]]); - // one more found + // If the current element belongs the the current column, + // add it in case it is also in the current row. + if ((positions[currentColumn] < outerPtr[currentColumn + 1]) + && (indexPtr[positions[currentColumn]] == currentRow)) { + addNextValue(currentRow, currentColumn, + valuePtr[positions[currentColumn]]); + // Remember that we found one more non-zero element. ++i; - // mark this position as "used" + // Mark this position as "used". ++positions[currentColumn]; } - // advance to next column + + // Now we can advance to the next column and also row, + // in case we just iterated through the last column. ++currentColumn; - if (currentColumn == eigen_col_count) { + if (currentColumn == colCount) { currentColumn = 0; ++currentRow; } @@ -227,335 +263,522 @@ class StaticSparseMatrix { } } - //! Linear Setter for matrix entry A_{row, col} to value /*! - Linear Setter function for matrix entry A_{row, col} to value. Must be called consecutively for each element in a row in ascending order of columns AND in ascending order of rows. - Diagonal entries may be set at any time. + * Sets the matrix element at the given row and column to the given value. + * NOTE: This is a linear setter. It must be called consecutively for each element, + * row by row *and* column by column. Only diagonal entries may be set at any time. + * @param row The row in which the matrix element is to be set. + * @param col The column in which the matrix element is to be set. + * @param value The value that is to be set. */ - void addNextValue(const uint_fast64_t row, const uint_fast64_t col, const T &value) { - if ((row > row_count) || (col > row_count) || (row == 0) || (col == 0)) { + void addNextValue(const uint_fast64_t row, const uint_fast64_t col, const T& value) { + // Check whether the given row and column positions are valid and throw + // error otherwise. + if ((row > rowCount) || (col > rowCount)) { triggerErrorState(); - pantheios::log_ERROR("StaticSparseMatrix::addNextValue: Throwing out_of_range: row or col not in 1 .. rows (is ", pantheios::integer(row), " x ", pantheios::integer(col), ", max is ", pantheios::integer(row_count), " x ", pantheios::integer(row_count), ")."); - throw mrmc::exceptions::out_of_range("StaticSparseMatrix::addNextValue: row or col not in 1 .. rows"); + pantheios::log_ERROR("StaticSparseMatrix::addNextValue: Throwing out_of_range: row or col not in 0 .. rows (is ", + pantheios::integer(row), " x ", pantheios::integer(col), ", max is ", + pantheios::integer(rowCount), " x ", pantheios::integer(rowCount), ")."); + throw mrmc::exceptions::out_of_range("StaticSparseMatrix::addNextValue: row or col not in 0 .. rows"); } - - if (row == col) { - diagonal_storage[row] = value; - } else { - if (row != last_row) { - for (uint_fast64_t i = last_row; i < row; ++i) { - row_indications[i] = current_size; + + if (row == col) { // Set a diagonal element. + diagonalStorage[row] = value; + } else { // Set a non-diagonal element. + // If we switched to another row, we have to adjust the missing + // entries in the row_indications array. + if (row != lastRow) { + for (uint_fast64_t i = lastRow + 1; i <= row; ++i) { + rowIndications[i] = currentSize; } - last_row = row; + lastRow = row; } - value_storage[current_size] = value; - column_indications[current_size] = col; + // Finally, set the element and increase the current size. + valueStorage[currentSize] = value; + columnIndications[currentSize] = col; - // Increment counter for stored elements - current_size++; + ++currentSize; } } + /* + * Finalizes the sparse matrix to indicate that initialization has been + * completed and the matrix may now be used. + */ void finalize() { + // Check whether it's safe to finalize the matrix and throw error + // otherwise. if (!isInitialized()) { triggerErrorState(); - pantheios::log_ERROR("StaticSparseMatrix::finalize: Throwing invalid state for internal state not Initialized (is ", pantheios::integer(internal_status)," - Already finalized?"); + pantheios::log_ERROR("StaticSparseMatrix::finalize: Throwing invalid state for internal state not Initialized (is ", + pantheios::integer(internalStatus), " - Already finalized?"); throw mrmc::exceptions::invalid_state("StaticSparseMatrix::finalize: Invalid state for internal state not Initialized - Already finalized?"); - } else if (current_size != non_zero_entry_count) { + } else if (currentSize != nonZeroEntryCount) { triggerErrorState(); pantheios::log_ERROR("StaticSparseMatrix::finalize: Throwing invalid_state: Wrong call count for addNextValue"); throw mrmc::exceptions::invalid_state("StaticSparseMatrix::finalize: Wrong call count for addNextValue"); } else { - if (last_row != row_count) { - for (uint_fast64_t i = last_row; i < row_count; ++i) { - row_indications[i] = current_size; + // Fill in the missing entries in the row_indications array. + // (Can happen because of empty rows at the end.) + if (lastRow != rowCount) { + for (uint_fast64_t i = lastRow + 1; i <= rowCount; ++i) { + rowIndications[i] = currentSize; } } - row_indications[row_count] = non_zero_entry_count; + // Set a sentinel element at the last position of the row_indications + // array. This eases iteration work, as now the indices of row i + // are always between row_indications[i] and row_indications[i + 1], + // also for the first and last row. + rowIndications[rowCount + 1] = nonZeroEntryCount; setState(MatrixStatus::ReadReady); } } - //! Getter for saving matrix entry A_{row,col} to target /*! - Getter function for the matrix. This function does not check the internal status for errors for performance reasons. - \param row 1-based index of the requested row - \param col 1-based index of the requested column - \param target pointer to where the result will be stored - \return True iff the value was set, false otherwise. On false, 0 will be written to *target. - */ + * Gets the matrix element at the given row and column to the given value. + * NOTE: This function does not check the internal status for errors for performance reasons. + * @param row The row in which the element is to be read. + * @param col The column in which the element is to be read. + * @param target A pointer to the memory location where the read content is + * to be put. + * @return True iff the value is set in the matrix, false otherwise. + * On false, 0 will be written to *target. + */ inline bool getValue(uint_fast64_t row, uint_fast64_t col, T* const target) { - - if (row == col) { - // storage is row_count + 1 large for direct access without the -1 - *target = diagonal_storage[row]; - return true; - } - - if ((row > row_count) || (col > row_count) || (row == 0) || (col == 0)) { - pantheios::log_ERROR("StaticSparseMatrix::getValue: row or col not in 1 .. rows (is ", pantheios::integer(row), " x ", pantheios::integer(col), ", max is ", pantheios::integer(row_count), " x ", pantheios::integer(row_count), ")."); - throw mrmc::exceptions::out_of_range("StaticSparseMatrix::getValue: row or col not in 1 .. rows"); + // Check for illegal access indices. + if ((row > rowCount) || (col > rowCount)) { + pantheios::log_ERROR("StaticSparseMatrix::getValue: row or col not in 0 .. rows (is ", pantheios::integer(row), " x ", + pantheios::integer(col), ", max is ", pantheios::integer(rowCount), " x ", pantheios::integer(rowCount), ")."); + throw mrmc::exceptions::out_of_range("StaticSparseMatrix::getValue: row or col not in 0 .. rows"); return false; } - uint_fast64_t row_start = row_indications[row - 1]; - uint_fast64_t row_end = row_indications[row]; + // Read elements on the diagonal directly. + if (row == col) { + *target = diagonalStorage[row]; + return true; + } - while (row_start < row_end) { - if (column_indications[row_start] == col) { - *target = value_storage[row_start]; + // In case the element is not on the diagonal, we have to iterate + // over the accessed row to find the element. + uint_fast64_t rowStart = rowIndications[row]; + uint_fast64_t rowEnd = rowIndications[row + 1]; + while (rowStart < rowEnd) { + // If the lement is found, write the content to the specified + // position and return true. + if (columnIndications[rowStart] == col) { + *target = valueStorage[rowStart]; return true; } - if (column_indications[row_start] > col) { + // If the column of the current element is already larger than the + // requested column, the requested element cannot be contained + // in the matrix and we may therefore stop searching. + if (columnIndications[rowStart] > col) { break; } - row_start++; + ++rowStart; } + // Set 0 as the content and return false in case the element was not found. *target = 0; return false; } + /*! + * Returns the number of rows of the matrix. + */ uint_fast64_t getRowCount() const { - return row_count; + return rowCount; } + /*! + * Returns a pointer to the value storage of the matrix. This storage does + * *not* include elements on the diagonal. + * @return A pointer to the value storage of the matrix. + */ T* getStoragePointer() const { - return value_storage; + return valueStorage; } + /*! + * Returns a pointer to the storage of elements on the diagonal. + * @return A pointer to the storage of elements on the diagonal. + */ T* getDiagonalStoragePointer() const { - return diagonal_storage; + return diagonalStorage; } + /*! + * Returns a pointer to the array that stores the start indices of non-zero + * entries in the value storage for each row. + * @return A pointer to the array that stores the start indices of non-zero + * entries in the value storage for each row. + */ uint_fast64_t* getRowIndicationsPointer() const { - return row_indications; + return rowIndications; } + /*! + * Returns a pointer to an array that stores the column of each non-zero + * element that is not on the diagonal. + * @return A pointer to an array that stores the column of each non-zero + * element that is not on the diagonal. + */ uint_fast64_t* getColumnIndicationsPointer() const { - return column_indications; + return columnIndications; } + /*! + * Checks whether the internal status of the matrix makes it ready for + * reading access. + * @return True iff the internal status of the matrix makes it ready for + * reading access. + */ bool isReadReady() { - return (internal_status == MatrixStatus::ReadReady); + return (internalStatus == MatrixStatus::ReadReady); } + /*! + * Checks whether the matrix was initialized previously. The matrix may + * still require to be finalized, even if this check returns true. + * @return True iff the matrix was initialized previously. + */ bool isInitialized() { - return (internal_status == MatrixStatus::Initialized || internal_status == MatrixStatus::ReadReady); + return (internalStatus == MatrixStatus::Initialized || internalStatus == MatrixStatus::ReadReady); } + /*! + * Returns the internal state of the matrix. + * @return The internal state of the matrix. + */ MatrixStatus getState() { - return internal_status; + return internalStatus; } + /*! + * Checks whether the internal state of the matrix signals an error. + * @return True iff the internal state of the matrix signals an error. + */ bool hasError() { - return (internal_status == MatrixStatus::Error); + return (internalStatus == MatrixStatus::Error); } - //! Converts this matrix to an equivalent sparse matrix in Eigens format. /*! - Exports this sparse matrix to Eigens SparseMatrix format. - Required this matrix to be in the ReadReady state. - - @return The Eigen SparseMatrix - + * Exports this sparse matrix to Eigens sparse matrix format. + * NOTE: this requires this matrix to be in the ReadReady state. + * @return The sparse matrix in Eigen format. */ Eigen::SparseMatrix* toEigenSparseMatrix() { - int_fast32_t eigenRows = static_cast(row_count); - Eigen::SparseMatrix* mat = new Eigen::SparseMatrix(eigenRows, eigenRows); - + // Check whether it is safe to export this matrix. if (!isReadReady()) { triggerErrorState(); - pantheios::log_ERROR("StaticSparseMatrix::toEigenSparseMatrix: Throwing invalid state for internal state not ReadReady (is ", pantheios::integer(internal_status),")."); + pantheios::log_ERROR("StaticSparseMatrix::toEigenSparseMatrix: Throwing invalid state for internal state not ReadReady (is ", pantheios::integer(internalStatus), ")."); throw mrmc::exceptions::invalid_state("StaticSparseMatrix::toEigenSparseMatrix: Invalid state for internal state not ReadReady."); } else { + // Create a + int_fast32_t eigenRows = static_cast(rowCount); + Eigen::SparseMatrix* mat = new Eigen::SparseMatrix(eigenRows, eigenRows); + + // There are two ways of converting this matrix to Eigen's format. + // 1. Compute a list of triplets (row, column, value) for all + // non-zero elements and pass it to Eigen to create a sparse matrix. + // 2. Tell Eigen to reserve the average number of non-zero elements + // per row in advance and insert the values via a call to Eigen's + // insert method then. As the given reservation number is only an + // estimate, the actual number may be different and Eigen may have + // to shift a lot. + // In most cases, the second alternative is faster (about 1/2 of the + // first, but if there are "heavy" rows that are several times larger + // than an average row, the other solution might be faster. + // The desired conversion method may be set by an appropriate define. # ifdef MRMC_USE_TRIPLETCONVERT + + // FIXME: Wouldn't it be more efficient to add the elements in + // order including the diagonal elements? Otherwise, Eigen has to + // perform some sorting. + + // Prepare the triplet storage. typedef Eigen::Triplet IntTriplet; std::vector tripletList; - tripletList.reserve(non_zero_entry_count + row_count); - - uint_fast64_t row_start; - uint_fast64_t row_end; - for (uint_fast64_t row = 1; row <= row_count; ++row) { - row_start = row_indications[row - 1]; - row_end = row_indications[row]; - while (row_start < row_end) { - tripletList.push_back(IntTriplet(row - 1, column_indications[row_start] - 1, value_storage[row_start])); - ++row_start; + tripletList.reserve(nonZeroEntryCount + rowCount); + + // First, iterate over all elements that are not on the diagonal + // and add the corresponding triplet. + uint_fast64_t rowStart; + uint_fast64_t rowEnd; + for (uint_fast64_t row = 0; row <= rowCount; ++row) { + rowStart = rowIndications[row]; + rowEnd = rowIndications[row + 1]; + while (rowStart < rowEnd) { + tripletList.push_back(IntTriplet(row, columnIndications[rowStart], valueStorage[rowStart])); + ++rowStart; } } - for (uint_fast64_t i = 1; i <= row_count; ++i) { - tripletList.push_back(IntTriplet(i - 1, i - 1, diagonal_storage[i])); + // Then add the elements on the diagonal. + for (uint_fast64_t i = 0; i <= rowCount; ++i) { + tripletList.push_back(IntTriplet(i, i, diagonalStorage[i])); } + // Let Eigen create a matrix from the given list of triplets. mat->setFromTriplets(tripletList.begin(), tripletList.end()); # else // NOT MRMC_USE_TRIPLETCONVERT - // In most cases, this is faster (about 1/2 of the above). But if there are "heavy" rows that are several times larger than an average row, the other solution might be faster. - - mat->reserve(Eigen::VectorXi::Constant(static_cast(mat->outerSize()), static_cast((non_zero_entry_count + row_count) / mat->outerSize()))); - - uint_fast64_t row_start; - uint_fast64_t row_end; - for (uint_fast64_t row = 1; row <= row_count; ++row) { - row_start = row_indications[row - 1]; - row_end = row_indications[row]; - - // insert the diagonal entry - mat->insert(row - 1, row - 1) = diagonal_storage[row]; - - while (row_start < row_end) { - //tripletList.push_back(IntTriplet(row - 1, column_indications[row_start] - 1, value_storage[row_start])); - mat->insert(row - 1, column_indications[row_start] - 1) = value_storage[row_start]; - ++row_start; + // Reserve the average number of non-zero elements per row for each + // row. + mat->reserve(Eigen::VectorXi::Constant(eigenRows, static_cast((nonZeroEntryCount + rowCount) / eigenRows))); + + // Iterate over the all non-zero elements in this matrix and add + // them to the matrix individually. + uint_fast64_t rowStart; + uint_fast64_t rowEnd; + 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]; + + // Insert the elements that are not on the diagonal + while (rowStart < rowEnd) { + mat->insert(row, columnIndications[rowStart]) = valueStorage[rowStart]; + ++rowStart; } } # endif // MRMC_USE_TRIPLETCONVERT + // Make the matrix compressed, i.e. remove remaining zero-entries. mat->makeCompressed(); + + return mat; } - - return mat; + + // This point can never be reached as both if-branches end in a return + // statement. + return nullptr; } - //! Returns the exact count of explicit entries in the sparse matrix. /*! - Retuns the exact count of explicit entries in the sparse matrix. While it is called "nonZero" count, the fields may of course be 0 (the 0-value of T). - @returns explicit entry count in the matrix + * 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. */ uint_fast64_t getNonZeroEntryCount() const { - return non_zero_entry_count; + return nonZeroEntryCount; } - //! Converts a state into a final state. /*! - This function allows for states to be made final. This means that all entries in row "state" will be changed to 0 but for the one on the diagonal, which will be set to 1 to create a loop on itself. - @param state The number of the state to be converted. Must be in 1 <= state <= rows. - @returns Whether the conversion was successful. + * This function makes the given state absorbing. This means that all + * entries in its row will be changed to 0 and the value 1 will be written + * to the element on the diagonal. + * @param state The state to be made absorbing. + * @returns True iff the operation was successful. */ - bool makeStateFinal(const uint_fast64_t state) { - if ((state > row_count) || (state == 0)) { - pantheios::log_ERROR("StaticSparseMatrix::makeStateFinal: state not in 1 .. rows (is ", pantheios::integer(state), ", max is ", pantheios::integer(row_count), ")."); - throw mrmc::exceptions::out_of_range("StaticSparseMatrix::makeStateFinal: state not in 1 .. rows"); + bool makeStateAbsorbing(const uint_fast64_t state) { + // Check whether the accessed state exists. + if (state > rowCount) { + pantheios::log_ERROR("StaticSparseMatrix::makeStateFinal: state not in 0 .. rows (is ", pantheios::integer(state), ", max is ", pantheios::integer(rowCount), ")."); + throw mrmc::exceptions::out_of_range("StaticSparseMatrix::makeStateFinal: state not in 0 .. rows"); return false; } - uint_fast64_t row_start = row_indications[state - 1]; - uint_fast64_t row_end = row_indications[state]; + // Iterate over the elements in the row that are not on the diagonal + // and set them to zero. + uint_fast64_t rowStart = rowIndications[state]; + uint_fast64_t rowEnd = rowIndications[state + 1]; - - while (row_start < row_end) { - value_storage[row_start] = mrmc::misc::constGetZero(value_storage); - row_start++; + while (rowStart < rowEnd) { + valueStorage[rowStart] = mrmc::misc::constGetZero(valueStorage); + ++rowStart; } - diagonal_storage[state] = mrmc::misc::constGetOne(diagonal_storage); + // Set the element on the diagonal to one. + diagonalStorage[state] = mrmc::misc::constGetOne(diagonalStorage); return true; } - /*! - * Returns the size of the matrix in memory measured in bytes. - * @return The size of the matrix in memory measured in bytes. - */ - uint_fast64_t getSizeInMemory() { - uint_fast64_t size = sizeof(*this); - size += sizeof(T) * non_zero_entry_count; // add value_storage size - size += sizeof(T) * (row_count + 1); // add diagonal_storage size - size += sizeof(uint_fast64_t) * non_zero_entry_count; // add column_indications size - size += sizeof(uint_fast64_t) * (row_count + 1); // add row_indications size - return size; - } - - private: - uint_fast64_t current_size; - - uint_fast64_t row_count; - uint_fast64_t non_zero_entry_count; - uint_fast64_t last_row; - - /*! Array containing all non-zero values, apart from the diagonal entries */ - T* value_storage; - /*! Array containing all diagonal values */ - T* diagonal_storage; - - /*! Array containing the column number of the corresponding value_storage entry */ - uint_fast64_t* column_indications; - /*! Array containing the row boundaries of valueStorage */ - uint_fast64_t* row_indications; - - /*! Internal status enum, 0 for constructed, 1 for initialized and 2 for finalized, -1 on errors */ - MatrixStatus internal_status; - - /*! Sets the internal status to Error */ + /*! + * Returns the size of the matrix in memory measured in bytes. + * @return The size of the matrix in memory measured in bytes. + */ + uint_fast64_t getSizeInMemory() { + uint_fast64_t size = sizeof(*this); + // Add value_storage size. + size += sizeof(T) * nonZeroEntryCount; + // Add diagonal_storage size. + size += sizeof(T) * (rowCount + 1); + // Add column_indications size. + size += sizeof(uint_fast64_t) * nonZeroEntryCount; + // Add row_indications size. + size += sizeof(uint_fast64_t) * (rowCount + 1); + return size; + } + + /*! + * Returns an iterator to the columns of the non-zero entries of the given + * row that are not on the diagonal. + * @param row The row whose columns the iterator will return. + * @return An iterator to the columns of the non-zero entries of the given + * row that are not on the diagonal. + */ + constIndexIterator beginConstColumnNoDiagIterator(uint_fast64_t row) const { + return this->columnIndications + this->rowIndications[row]; + } + + /*! + * Returns an iterator referring to the element after the given row. + * @param row The row for which the iterator should point to the past-the-end + * element. + */ + constIndexIterator endConstColumnNoDiagIterator(uint_fast64_t row) const { + return this->columnIndications + this->rowIndications[row + 1]; + } + +private: + + /*! + * The number of rows of the matrix. + */ + uint_fast64_t rowCount; + + /*! + * The number of non-zero elements that are not on the diagonal. + */ + uint_fast64_t nonZeroEntryCount; + + /*! + * Stores all non-zero values that are not on the diagonal. + */ + T* valueStorage; + + /*! + * Stores all elements on the diagonal, even the ones that are zero. + */ + T* diagonalStorage; + + /*! + * Stores the column for each non-zero element that is not on the diagonal. + */ + uint_fast64_t* columnIndications; + + /*! + * Array containing the boundaries (indices) in the value_storage array + * for each row. All elements of value_storage with indices between the + * i-th and the (i+1)-st element of this array belong to row i. + */ + uint_fast64_t* rowIndications; + + /*! + * The internal status of the matrix. + */ + MatrixStatus internalStatus; + + /*! + * Stores the current number of non-zero elements that have been added to + * the matrix. Used for correctly inserting elements in the matrix. + */ + uint_fast64_t currentSize; + + /*! + * Stores the row in which the last element was inserted. Used for correctly + * inserting elements in the matrix . + */ + uint_fast64_t lastRow; + + /*! + * Sets the internal status to signal an error. + */ void triggerErrorState() { setState(MatrixStatus::Error); } /*! - Sets the internal status to new_state iff the current state is not the Error state. - @param new_state the new state to be switched to - */ + * Sets the internal status to the given state if the current state is not + * the error state. + * @param new_state The new state to be switched to. + */ void setState(const MatrixStatus new_state) { - internal_status = (internal_status == MatrixStatus::Error) ? internal_status : new_state; + internalStatus = (internalStatus == MatrixStatus::Error) ? internalStatus : new_state; } /*! - Prepares the internal CSR storage. - Requires non_zero_entry_count and row_count to be set. - @return true on success, false otherwise (allocation failed). + * Prepares the internal CSR storage. For this, it requires + * non_zero_entry_count and row_count to be set correctly. + * @return True on success, false otherwise (allocation failed). */ bool prepareInternalStorage() { - value_storage = new (std::nothrow) T[non_zero_entry_count](); + // Set up the arrays for the elements that are not on the diagonal. + valueStorage = new (std::nothrow) T[nonZeroEntryCount](); + columnIndications = new (std::nothrow) uint_fast64_t[nonZeroEntryCount](); - column_indications = new (std::nothrow) uint_fast64_t[non_zero_entry_count](); + // Set up the row_indications array and reserve one element more than + // there are rows in order to put a sentinel element at the end, + // which eases iteration process. + rowIndications = new (std::nothrow) uint_fast64_t[rowCount + 1](); - row_indications = new (std::nothrow) uint_fast64_t[row_count + 1](); + // Set up the array for the elements on the diagonal. + diagonalStorage = new (std::nothrow) T[rowCount](); - // row_count + 1 so that access with 1-based indices can be direct without the overhead of a -1 each time - diagonal_storage = new (std::nothrow) T[row_count + 1](); - - return ((value_storage != NULL) && (column_indications != NULL) && (row_indications != NULL) && (diagonal_storage != NULL)); + // Return whether all the allocations could be made without error. + return ((valueStorage != NULL) && (columnIndications != NULL) + && (rowIndications != NULL) && (diagonalStorage != NULL)); } - //! - - - template + /*! + * Helper function to determine whether the given Eigen matrix is in RowMajor + * format. Always returns true, but is overloaded, so the compiler will + * only call it in case the Eigen matrix is in RowMajor format. + * @return True. + */ + template bool isEigenRowMajor(Eigen::SparseMatrix<_Scalar, Eigen::RowMajor, _Index>) { return true; } - - template - bool isEigenRowMajor(Eigen::SparseMatrix<_Scalar, Eigen::ColMajor, _Index>) { + + /*! + * Helper function to determine whether the given Eigen matrix is in RowMajor + * format. Always returns false, but is overloaded, so the compiler will + * only call it in case the Eigen matrix is in ColMajor format. + * @return False. + */ + template + bool isEigenRowMajor( + Eigen::SparseMatrix<_Scalar, Eigen::ColMajor, _Index>) { return false; } + /*! + * Helper function to determine the number of non-zero elements that are + * not on the diagonal of the given Eigen matrix. + * @param eigen_sparse_matrix The Eigen matrix to analyze. + * @return The number of non-zero elements that are not on the diagonal of + * the given Eigen matrix. + */ template - _Index getEigenSparseMatrixCorrectNonZeroEntryCount(const Eigen::SparseMatrix<_Scalar, _Options, _Index> &eigen_sparse_matrix) { + _Index getEigenSparseMatrixCorrectNonZeroEntryCount(const Eigen::SparseMatrix<_Scalar, _Options, _Index>& eigen_sparse_matrix) { const _Index* indexPtr = eigen_sparse_matrix.innerIndexPtr(); const _Index* outerPtr = eigen_sparse_matrix.outerIndexPtr(); const _Index entryCount = eigen_sparse_matrix.nonZeros(); const _Index outerCount = eigen_sparse_matrix.outerSize(); - - uint_fast64_t diag_non_zeros = 0; - - // for RowMajor, row is the current Row and col the column - // for ColMajor, row is the current Col and col the row + + uint_fast64_t diagNonZeros = 0; + + // For RowMajor, row is the current row and col the column and for + // ColMajor, row is the current column and col the row, but this is + // not important as we are only looking for elements on the diagonal. _Index innerStart = 0; - _Index innerEnd = 0; - _Index innerMid = 0; + _Index innerEnd = 0; + _Index innerMid = 0; for (_Index row = 0; row < outerCount; ++row) { innerStart = outerPtr[row]; - innerEnd = outerPtr[row + 1] - 1; - - // Now with super fancy binary search, deferred equality detection + innerEnd = outerPtr[row + 1] - 1; + + // Now use binary search (but defer equality detection). while (innerStart < innerEnd) { innerMid = innerStart + ((innerEnd - innerStart) / 2); @@ -566,13 +789,13 @@ class StaticSparseMatrix { } } + // Check whether we have found an element on the diagonal. if ((innerStart == innerEnd) && (indexPtr[innerStart] == row)) { - // found a diagonal entry - ++diag_non_zeros; + ++diagNonZeros; } } - return static_cast<_Index>(entryCount - diag_non_zeros); + return static_cast<_Index>(entryCount - diagNonZeros); } }; diff --git a/test/sparse/static_sparse_matrix_test.cpp b/test/sparse/static_sparse_matrix_test.cpp index 167f221f5..552534a23 100644 --- a/test/sparse/static_sparse_matrix_test.cpp +++ b/test/sparse/static_sparse_matrix_test.cpp @@ -29,10 +29,10 @@ TEST(StaticSparseMatrixTest, addNextValueTest) { ASSERT_NO_THROW(ssm->initialize(1)); ASSERT_EQ(ssm->getState(), mrmc::sparse::StaticSparseMatrix::MatrixStatus::Initialized); - ASSERT_THROW(ssm->addNextValue(0, 1, 1), mrmc::exceptions::out_of_range); + ASSERT_THROW(ssm->addNextValue(-1, 1, 1), mrmc::exceptions::out_of_range); ASSERT_EQ(ssm->getState(), mrmc::sparse::StaticSparseMatrix::MatrixStatus::Error); - ASSERT_THROW(ssm->addNextValue(1, 0, 1), mrmc::exceptions::out_of_range); + ASSERT_THROW(ssm->addNextValue(1, -1, 1), mrmc::exceptions::out_of_range); ASSERT_EQ(ssm->getState(), mrmc::sparse::StaticSparseMatrix::MatrixStatus::Error); ASSERT_THROW(ssm->addNextValue(6, 1, 1), mrmc::exceptions::out_of_range); @@ -78,7 +78,7 @@ TEST(StaticSparseMatrixTest, Test) { int position_row[50] = { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, - 2, 2, 2, 2, /* first row empty, one full row � 25 minus the diagonal entry */ + 2, 2, 2, 2, /* first row empty, one full row ��� 25 minus the diagonal entry */ 4, 4, /* one empty row, then first and last column */ 13, 13, 13, 13, /* a few empty rows, middle columns */ 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, @@ -151,10 +151,10 @@ TEST(StaticSparseMatrixTest, ConversionFromDenseEigen_ColMajor_SparseMatrixTest) ASSERT_EQ(ssm->getState(), mrmc::sparse::StaticSparseMatrix::MatrixStatus::ReadReady); int target = -1; - for (int row = 1; row <= 10; ++row) { - for (int col = 1; col <= 10; ++col) { + for (int row = 0; row < 10; ++row) { + for (int col = 0; col < 10; ++col) { ASSERT_TRUE(ssm->getValue(row, col, &target)); - ASSERT_EQ(target, (row - 1) * 10 + col - 1); + ASSERT_EQ(target, row * 10 + col); } } } @@ -181,10 +181,10 @@ TEST(StaticSparseMatrixTest, ConversionFromDenseEigen_RowMajor_SparseMatrixTest) ASSERT_EQ(ssm->getState(), mrmc::sparse::StaticSparseMatrix::MatrixStatus::ReadReady); int target = -1; - for (int row = 1; row <= 10; ++row) { - for (int col = 1; col <= 10; ++col) { + for (int row = 0; row < 10; ++row) { + for (int col = 0; col < 10; ++col) { ASSERT_TRUE(ssm->getValue(row, col, &target)); - ASSERT_EQ(target, (row - 1) * 10 + col - 1); + ASSERT_EQ(target, row * 10 + col); } } } @@ -231,7 +231,7 @@ TEST(StaticSparseMatrixTest, ConversionFromSparseEigen_ColMajor_SparseMatrixTest int target = -1; for (auto &coeff: tripletList) { - ASSERT_TRUE(ssm->getValue(coeff.row() + 1, coeff.col() + 1, &target)); + ASSERT_TRUE(ssm->getValue(coeff.row(), coeff.col(), &target)); ASSERT_EQ(target, coeff.value()); } } @@ -278,7 +278,7 @@ TEST(StaticSparseMatrixTest, ConversionFromSparseEigen_RowMajor_SparseMatrixTest int target = -1; for (auto &coeff: tripletList) { - ASSERT_TRUE(ssm->getValue(coeff.row() + 1, coeff.col() + 1, &target)); + ASSERT_TRUE(ssm->getValue(coeff.row(), coeff.col(), &target)); ASSERT_EQ(target, coeff.value()); } } @@ -294,9 +294,9 @@ TEST(StaticSparseMatrixTest, ConversionToSparseEigen_RowMajor_SparseMatrixTest) ASSERT_NO_THROW(ssm->initialize(100 - 10)); ASSERT_EQ(ssm->getState(), mrmc::sparse::StaticSparseMatrix::MatrixStatus::Initialized); - for (uint_fast32_t row = 1; row <= 10; ++row) { - for (uint_fast32_t col = 1; col <= 10; ++col) { - ASSERT_NO_THROW(ssm->addNextValue(row, col, values[(row - 1) * 10 + (col - 1)])); + for (uint_fast32_t row = 0; row < 10; ++row) { + for (uint_fast32_t col = 0; col < 10; ++col) { + ASSERT_NO_THROW(ssm->addNextValue(row, col, values[row * 10 + col])); } } ASSERT_EQ(ssm->getState(), mrmc::sparse::StaticSparseMatrix::MatrixStatus::Initialized); @@ -311,4 +311,4 @@ TEST(StaticSparseMatrixTest, ConversionToSparseEigen_RowMajor_SparseMatrixTest) ASSERT_EQ(values[row * 10 + col], esm->coeff(row, col)); } } -} \ No newline at end of file +}