diff --git a/src/dtmc/atomic_proposition.h b/src/dtmc/atomic_proposition.h index 001a92a4b..5beb41227 100644 --- a/src/dtmc/atomic_proposition.h +++ b/src/dtmc/atomic_proposition.h @@ -8,10 +8,7 @@ #include #include -#include "src/exceptions/invalid_state.h" -#include "src/exceptions/invalid_argument.h" -#include "src/exceptions/out_of_range.h" - +#include "src/vector/bitvector.h" namespace mrmc { @@ -29,47 +26,25 @@ class AtomicProposition { /*! \param nodeCount Amount of nodes that the DTMC has to label */ - AtomicProposition(uint_fast32_t nodeCount) { - node_array = NULL; - if (node_array == NULL) { - node_count = nodeCount; - int n = (int) std::ceil(nodeCount / 64.0); - node_array = new uint_fast64_t[n]; - //Initialization with 0 is crucial! - memset(node_array, 0, n*sizeof(uint_fast64_t)); - } + AtomicProposition(uint_fast32_t nodeCount) : nodes(nodeCount) { + // } ~AtomicProposition() { - if (node_array != NULL) { - delete[] node_array; - } + // } bool hasNodeLabel(uint_fast32_t nodeId) { - int bucket = static_cast(std::floor(nodeId / 64)); - int bucket_place = nodeId % 64; - // Taking the step with mask is crucial as we NEED a 64bit shift, not a 32bit one. - uint_fast64_t mask = 1; - mask = mask << bucket_place; - return ((mask & node_array[bucket]) == mask); + return nodes.get(nodeId); } void addLabelToNode(uint_fast32_t nodeId) { - int bucket = static_cast(std::floor(nodeId / 64)); - // Taking the step with mask is crucial as we NEED a 64bit shift, not a 32bit one. - // MSVC: C4334, use 1i64 or cast to __int64. - // result of 32-bit shift implicitly converted to 64 bits (was 64-bit shift intended?) - uint_fast64_t mask = 1; - mask = mask << (nodeId % 64); - node_array[bucket] |= mask; + nodes.set(nodeId, true); } private: - uint_fast32_t node_count; - - /*! Array containing the boolean bits for each node, 64bits/64nodes per element */ - uint_fast64_t* node_array; + /*! BitVector containing the boolean bits for each node */ + mrmc::vector::BitVector nodes; }; } // namespace dtmc diff --git a/src/mrmc-cpp.cpp b/src/mrmc-cpp.cpp index 666723a18..cbbeabb67 100644 --- a/src/mrmc-cpp.cpp +++ b/src/mrmc-cpp.cpp @@ -15,7 +15,8 @@ #include #include -#include +/* PlatformSTL Header Files */ +#include #include #include @@ -26,60 +27,88 @@ PANTHEIOS_EXTERN_C PAN_CHAR_T const PANTHEIOS_FE_PROCESS_IDENTITY[] = "mrmc-cpp" #include "src/sparse/static_sparse_matrix.h" #include "src/dtmc/atomic_proposition.h" #include "src/parser/read_lab_file.h" +#include "src/parser/read_tra_file.h" +#include "Eigen/Sparse" int main(int argc, char* argv[]) { // Logging init pantheios_be_file_setFilePath("log.all"); pantheios::log_INFORMATIONAL("MRMC-Cpp started."); - mrmc::dtmc::labeling *lab; + if (argc < 2) { + 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"); - time_t start = std::clock(); + counter.stop(); - lab = mrmc::parser::read_lab_file(20302, "csl_unbounded_until_sim_06.lab"); + 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; - time_t end = std::clock(); - printf("Time needed was %f", ((float)end-start)/CLOCKS_PER_SEC); - delete lab; + counter.start(); + auto esm = ssm->toEigenSparseMatrix(); + counter.stop(); - /* - std::cout << "Hello, World." << std::endl; - std::cout << "This is MRMC-Cpp Version " << MRMC_CPP_VERSION_MAJOR << "." << MRMC_CPP_VERSION_MINOR << std::endl; + 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 *ssm = new mrmc::sparse::StaticSparseMatrix(10, 10); + ssmBack = new mrmc::sparse::StaticSparseMatrix(esm->rows()); - ssm->initialize(); - ssm->addNextValue(2, 3, 1); - ssm->addNextValue(2, 4, 2); - ssm->addNextValue(2, 5, 3); - ssm->addNextValue(2, 6, 4); - ssm->addNextValue(2, 7, 5); + counter.start(); + ssmBack->initialize(*esm); + counter.stop(); - ssm->addNextValue(3, 4, 6); - ssm->addNextValue(3, 5, 7); - ssm->addNextValue(3, 6, 8); - ssm->addNextValue(3, 7, 9); - ssm->addNextValue(3, 8, 10); + 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; - ssm->finalize(); + delete ssm; + delete ssmBack; + + std::cout << std::endl; + std::cout << ":: C-style vs. C++ Style Array init Showdown ::" << std::endl; - int target; + uint_fast64_t *iArray = NULL; + uint_fast32_t i = 0; + uint_fast32_t j = 0; - for (int row = 1; row <= 10; ++row) { - for (int col = 1; col <= 10; ++col) { - if (ssm->getValue(row, col, &target)) { - printf(" T%i", target); - } else { - printf(" _%i", target); - } + 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; } - printf("\n"); + free(iArray); } + counter.stop(); - delete ssm; - */ + 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/sparse/static_sparse_matrix.h b/src/sparse/static_sparse_matrix.h index 7a2b2ee52..bc5abcb99 100644 --- a/src/sparse/static_sparse_matrix.h +++ b/src/sparse/static_sparse_matrix.h @@ -13,7 +13,6 @@ #include "src/exceptions/out_of_range.h" #include "Eigen/Sparse" -#include "src/sparse/eigen_sparse_additions.h" namespace mrmc { @@ -49,7 +48,7 @@ class StaticSparseMatrix { /*! \param rows Row-Count and therefore column-count of the square matrix */ - StaticSparseMatrix(uint_fast32_t rows) { + 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; @@ -66,6 +65,37 @@ class StaticSparseMatrix { //initialize(rows, non_zero_entries); } + //! 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 + */ + 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."); + if (!ssm.hasError()) { + pantheios::log_ERROR("StaticSparseMatrix::CopyCtor: Throwing invalid_argument: Can not Copy from matrix in Error state."); + throw mrmc::exceptions::invalid_argument(); + } else { + if (!prepareInternalStorage()) { + pantheios::log_ERROR("StaticSparseMatrix::CopyCtor: Throwing bad_alloc: memory allocation failed."); + throw std::bad_alloc(); + } else { + for (uint_fast64_t i = 0; i < non_zero_entry_count; ++i) { + // use T() to force use of the copy constructor for complex T types + value_storage[i] = T(ssm.value_storage[i]); + } + for (uint_fast64_t i = 0; i <= row_count; ++i) { + // same here + diagonal_storage[i] = T(ssm.diagonal_storage[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)); + } + } + } + ~StaticSparseMatrix() { setState(MatrixStatus::UnInitialized); if (value_storage != NULL) { @@ -93,7 +123,7 @@ class StaticSparseMatrix { 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) */ - void initialize(uint_fast32_t non_zero_entries) { + void initialize(uint_fast64_t non_zero_entries) { if (internal_status != MatrixStatus::UnInitialized) { triggerErrorState(); pantheios::log_ERROR("StaticSparseMatrix::initialize: Throwing invalid state for status flag != 0 (is ", pantheios::integer(internal_status)," - Already initialized?"); @@ -126,8 +156,8 @@ class StaticSparseMatrix { 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! */ - template - void initialize(const Eigen::SparseMatrix<_Scalar, _Options, _Index> &eigen_sparse_matrix) { + template + void initialize(const Eigen::SparseMatrix &eigen_sparse_matrix) { if (!eigen_sparse_matrix.isCompressed()) { triggerErrorState(); pantheios::log_ERROR("StaticSparseMatrix::initialize: Throwing invalid_argument: eigen_sparse_matrix is not in Compressed form."); @@ -145,28 +175,28 @@ class StaticSparseMatrix { // easy case, we can simply copy the data // RowMajor: Easy, ColMajor: Hmm. But how to detect? const T* valuePtr = eigen_sparse_matrix.valuePtr(); - const int_fast32_t* indexPtr = eigen_sparse_matrix.innerIndexPtr(); - const int_fast32_t* outerPtr = eigen_sparse_matrix.outerIndexPtr(); + const _Index* indexPtr = eigen_sparse_matrix.innerIndexPtr(); + const _Index* outerPtr = eigen_sparse_matrix.outerIndexPtr(); - const int_fast32_t entryCount = eigen_sparse_matrix.nonZeros(); - const int_fast32_t outerCount = eigen_sparse_matrix.outerSize(); + 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 (int row = 1; row <= outerCount; ++row) { - for (int col = outerPtr[row - 1]; col < outerPtr[row]; ++col) { + for (_Index row = 1; row <= outerCount; ++row) { + for (_Index col = outerPtr[row - 1]; col < outerPtr[row]; ++col) { addNextValue(row, indexPtr[col] + 1, valuePtr[col]); } } } else { // temp copies, anyone? - const int eigen_col_count = eigen_sparse_matrix.cols(); - const int eigen_row_count = eigen_sparse_matrix.rows(); + 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 - int_fast32_t* positions = new int_fast32_t[eigen_col_count](); - for (int i = 0; i < eigen_col_count; ++i) { + _Index* positions = new _Index[eigen_col_count](); + for (_Index i = 0; i < eigen_col_count; ++i) { positions[i] = outerPtr[i]; } @@ -188,6 +218,7 @@ class StaticSparseMatrix { ++currentRow; } } + delete[] positions; } setState(MatrixStatus::Initialized); } @@ -198,7 +229,7 @@ class StaticSparseMatrix { 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. */ - void addNextValue(const uint_fast32_t row, const uint_fast32_t col, const T &value) { + 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)) { 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), ")."); @@ -209,7 +240,7 @@ class StaticSparseMatrix { diagonal_storage[row] = value; } else { if (row != last_row) { - for (uint_fast32_t i = last_row; i < row; ++i) { + for (uint_fast64_t i = last_row; i < row; ++i) { row_indications[i] = current_size; } last_row = row; @@ -234,7 +265,7 @@ class StaticSparseMatrix { throw mrmc::exceptions::invalid_state("StaticSparseMatrix::finalize: Wrong call count for addNextValue"); } else { if (last_row != row_count) { - for (uint_fast32_t i = last_row; i < row_count; ++i) { + for (uint_fast64_t i = last_row; i < row_count; ++i) { row_indications[i] = current_size; } } @@ -253,7 +284,7 @@ class StaticSparseMatrix { \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. */ - inline bool getValue(uint_fast32_t row, uint_fast32_t col, T* const 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 @@ -267,8 +298,8 @@ class StaticSparseMatrix { return false; } - uint_fast32_t row_start = row_indications[row - 1]; - uint_fast32_t row_end = row_indications[row]; + uint_fast64_t row_start = row_indications[row - 1]; + uint_fast64_t row_end = row_indications[row]; while (row_start < row_end) { if (column_indications[row_start] == col) { @@ -285,7 +316,7 @@ class StaticSparseMatrix { return false; } - uint_fast32_t getRowCount() const { + uint_fast64_t getRowCount() const { return row_count; } @@ -317,21 +348,24 @@ class StaticSparseMatrix { @return The Eigen SparseMatrix */ - Eigen::SparseMatrix toEigenSparseMatrix() { - Eigen::SparseMatrix mat(row_count, row_count); + Eigen::SparseMatrix* toEigenSparseMatrix() { + int_fast32_t eigenRows = static_cast(row_count); + Eigen::SparseMatrix* mat = new Eigen::SparseMatrix(eigenRows, eigenRows); if (!isReadReady()) { triggerErrorState(); pantheios::log_ERROR("StaticSparseMatrix::toEigenSparseMatrix: Throwing invalid state for internal state not ReadReady (is ", pantheios::integer(internal_status),")."); throw mrmc::exceptions::invalid_state("StaticSparseMatrix::toEigenSparseMatrix: Invalid state for internal state not ReadReady."); } else { - typedef Eigen::Triplet IntTriplet; + +# ifdef MRMC_USE_TRIPLETCONVERT + typedef Eigen::Triplet IntTriplet; std::vector tripletList; tripletList.reserve(non_zero_entry_count + row_count); - uint_fast32_t row_start; - uint_fast32_t row_end; - for (uint_fast32_t row = 1; row <= row_count; ++row) { + 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) { @@ -339,88 +373,78 @@ class StaticSparseMatrix { ++row_start; } } - for (uint_fast32_t i = 1; i <= row_count; ++i) { - tripletList.push_back(IntTriplet(i, i, diagonal_storage[i])); + + for (uint_fast64_t i = 1; i <= row_count; ++i) { + tripletList.push_back(IntTriplet(i - 1, i - 1, diagonal_storage[i])); } - mat.setFromTriplets(tripletList.begin(), tripletList.end()); - mat.makeCompressed(); + 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]; + 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; + } + } +# endif // MRMC_USE_TRIPLETCONVERT + + mat->makeCompressed(); } return mat; } - //! Alternative way to initialize this matrix instead of using initialize(), addNextValue() and finalize(). + //! Returns the exact count of explicit entries in the sparse matrix. /*! - Initializes the matrix from the given eigen_sparse_matrix. Replaces the calls to initialize(), addNextValue() and finalize(). - Requires eigen_sparse_matrix to have at most a size of rows x rows and not more than non_zero_entries on non-diagonal fields. - To calculate the non-zero-entry count on only non-diagonal fields, you may use getEigenSparseMatrixCorrectNonZeroEntryCount(). - In most cases, it may be easier to use the alternative constructor StaticSparseMatrix(const Eigen::SparseMatrix eigen_sparse_matrix). - @param eigen_sparse_matrix the Eigen Sparse Matrix from which this matrix should be initalized. - @see getEigenSparseMatrixCorrectNonZeroEntryCount() - @see StaticSparseMatrix(const Eigen::SparseMatrix) + 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 */ - bool fromEigenSparseMatrix(const Eigen::SparseMatrix eigen_sparse_matrix) { - if (getState() != MatrixStatus::UnInitialized) { - triggerErrorState(); - pantheios::log_ERROR("StaticSparseMatrix::fromEigenSparseMatrix: Throwing invalid state for internal state not UnInitialized (is ", pantheios::integer(internal_status),")."); - throw mrmc::exceptions::invalid_state("StaticSparseMatrix::fromEigenSparseMatrix: Invalid state for internal state not UnInitialized."); - return false; - } - - int_fast32_t eigen_row_count = eigen_sparse_matrix.rows(); - int_fast32_t eigen_col_count = eigen_sparse_matrix.cols(); - - if ((eigen_row_count > row_count) || (eigen_col_count > row_count)) { - triggerErrorState(); - pantheios::log_ERROR("StaticSparseMatrix::fromEigenSparseMatrix: Throwing invalid argument for eigenSparseMatrix is too big to fit (is ", pantheios::integer(eigen_row_count)," x ", pantheios::integer(eigen_col_count), ", max is ", pantheios::integer(row_count)," x ", pantheios::integer(row_count), ")."); - throw mrmc::exceptions::invalid_argument("StaticSparseMatrix::fromEigenSparseMatrix: Invalid argument for eigenSparseMatrix is too big to fit."); - return false; - } - - uint_fast32_t eigen_non_zero_entries = mrmc::sparse::getEigenSparseMatrixCorrectNonZeroEntryCount(eigen_sparse_matrix); + uint_fast64_t getNonZeroEntryCount() const { + return non_zero_entry_count; + } - if (eigen_non_zero_entries > non_zero_entry_count) { - triggerErrorState(); - pantheios::log_ERROR("StaticSparseMatrix::fromEigenSparseMatrix: Throwing invalid argument for eigenSparseMatrix has too many non-zero entries to fit (is ", pantheios::integer(eigen_non_zero_entries),", max is ", pantheios::integer(non_zero_entry_count),")."); - throw mrmc::exceptions::invalid_argument("StaticSparseMatrix::fromEigenSparseMatrix: Invalid argument for eigenSparseMatrix has too many non-zero entries to fit."); + //! 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. + */ + 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"); return false; } - // make compressed - eigen_sparse_matrix.makeCompressed(); + uint_fast64_t row_start = row_indications[state - 1]; + uint_fast64_t row_end = row_indications[state]; - if (eigen_sparse_matrix.IsRowMajor()) { - // inner Index - int_fast32_t* eigenInnerIndex = eigen_sparse_matrix.innerIndexPtr(); - T* eigenValuePtr = eigen_sparse_matrix.valuePtr(); + while (row_start < row_end) { + value_storage[row_start] = getConstZero(); + row_start++; } - - - for (int k = 0; k < tempESM.outerSize(); ++k) { - for (SparseMatrix::InnerIterator it(tempESM, k); it; ++it) { - if (eigen_non_zero_entries >= non_zero_entry_count) { - // too many non zero entries for us. - } - addNextValue(it.row() - 1, it.col() - 1, it.value()); - if (it.row() != it.col()) { - ++eigen_non_zero_entries; - } - /*it.value(); - it.row(); // row index - it.col(); // col index (here it is equal to k) - it.index(); // inner index, here it is equal to it.row()*/ - } - } + diagonal_storage[state] = getConstOne(); + return true; } private: - uint_fast32_t current_size; + uint_fast64_t current_size; - uint_fast32_t row_count; - uint_fast32_t non_zero_entry_count; - uint_fast32_t last_row; + 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; @@ -428,9 +452,9 @@ class StaticSparseMatrix { T* diagonal_storage; /*! Array containing the column number of the corresponding value_storage entry */ - uint_fast32_t* column_indications; + uint_fast64_t* column_indications; /*! Array containing the row boundaries of valueStorage */ - uint_fast32_t* row_indications; + uint_fast64_t* row_indications; /*! Internal status enum, 0 for constructed, 1 for initialized and 2 for finalized, -1 on errors */ MatrixStatus internal_status; @@ -456,9 +480,9 @@ class StaticSparseMatrix { bool prepareInternalStorage() { value_storage = new (std::nothrow) T[non_zero_entry_count](); - column_indications = new (std::nothrow) uint_fast32_t[non_zero_entry_count](); + column_indications = new (std::nothrow) uint_fast64_t[non_zero_entry_count](); - row_indications = new (std::nothrow) uint_fast32_t[row_count + 1](); + row_indications = new (std::nothrow) uint_fast64_t[row_count + 1](); // 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](); @@ -466,6 +490,38 @@ class StaticSparseMatrix { return ((value_storage != NULL) && (column_indications != NULL) && (row_indications != NULL) && (diagonal_storage != NULL)); } + //! + + template + _Scalar constGetZero() const { + return _Scalar(0); + } + + template <> + int_fast32_t constGetZero() const { + return 0; + } + + template <> + double constGetZero() const { + return 0.0; + } + + template + _Scalar constGetOne() const { + return _Scalar(1); + } + + template<> + int_fast32_t constGetOne() const { + return 1; + } + + template<> + double constGetOne() const { + return 1.0; + } + template bool isEigenRowMajor(Eigen::SparseMatrix<_Scalar, Eigen::RowMajor, _Index>) { return true; @@ -477,21 +533,21 @@ class StaticSparseMatrix { } template - uint_fast32_t getEigenSparseMatrixCorrectNonZeroEntryCount(const Eigen::SparseMatrix<_Scalar, _Options, _Index> &eigen_sparse_matrix) { - const int_fast32_t* indexPtr = eigen_sparse_matrix.innerIndexPtr(); - const int_fast32_t* outerPtr = eigen_sparse_matrix.outerIndexPtr(); + _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 int_fast32_t entryCount = eigen_sparse_matrix.nonZeros(); - const int_fast32_t outerCount = eigen_sparse_matrix.outerSize(); + const _Index entryCount = eigen_sparse_matrix.nonZeros(); + const _Index outerCount = eigen_sparse_matrix.outerSize(); - uint_fast32_t diag_non_zeros = 0; + 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 - int_fast32_t innerStart = 0; - int_fast32_t innerEnd = 0; - int_fast32_t innerMid = 0; - for (int row = 0; row < outerCount; ++row) { + _Index innerStart = 0; + _Index innerEnd = 0; + _Index innerMid = 0; + for (_Index row = 0; row < outerCount; ++row) { innerStart = outerPtr[row]; innerEnd = outerPtr[row + 1] - 1; @@ -512,7 +568,7 @@ class StaticSparseMatrix { } } - return static_cast(entryCount - diag_non_zeros); + return static_cast<_Index>(entryCount - diag_non_zeros); } }; diff --git a/src/vector/bitvector.h b/src/vector/bitvector.h new file mode 100644 index 000000000..ec50d2870 --- /dev/null +++ b/src/vector/bitvector.h @@ -0,0 +1,152 @@ +#ifndef MRMC_VECTOR_BITVECTOR_H_ +#define MRMC_VECTOR_BITVECTOR_H_ + +#include +#include +#include +#include "boost/integer/integer_mask.hpp" + +#include +#include + +#include "src/exceptions/invalid_state.h" +#include "src/exceptions/invalid_argument.h" +#include "src/exceptions/out_of_range.h" + +namespace mrmc { + +namespace vector { + +//! A Boolean Array +/*! + A bit vector for boolean fields or quick selection schemas on Matrix entries. + Does NOT perform index bound checks! + */ +class BitVector { + public: + //! Constructor + /*! + \param initial_length The initial size of the boolean Array. Can be changed later on via @resize() + */ + BitVector(uint_fast64_t initial_length) { + bucket_count = initial_length / 64; + if (initial_length % 64 != 0) { + ++bucket_count; + } + bucket_array = new uint_fast64_t[bucket_count](); + + // init all 0 + for (uint_fast64_t i = 0; i < bucket_count; ++i) { + bucket_array[i] = 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 + */ + BitVector(const BitVector &bv) : bucket_count(bv.bucket_count) + { + pantheios::log_DEBUG("BitVector::CopyCTor: Using Copy() Ctor."); + bucket_array = new uint_fast64_t[bucket_count](); + memcpy(bucket_array, bv.bucket_array, sizeof(uint_fast64_t) * bucket_count); + } + + ~BitVector() { + if (bucket_array != NULL) { + delete[] bucket_array; + } + } + + void resize(uint_fast64_t new_length) { + uint_fast64_t* tempArray = new uint_fast64_t[new_length](); + + // 64 bit/entries per uint_fast64_t + uint_fast64_t copySize = (new_length <= (bucket_count * 64)) ? (new_length/64) : (bucket_count); + memcpy(tempArray, bucket_array, sizeof(uint_fast64_t) * copySize); + + bucket_count = new_length / 64; + if (new_length % 64 != 0) { + ++bucket_count; + } + + delete[] bucket_array; + bucket_array = tempArray; + } + + void set(const uint_fast64_t index, const bool value) { + uint_fast64_t bucket = index / 64; + // Taking the step with mask is crucial as we NEED a 64bit shift, not a 32bit one. + // MSVC: C4334, use 1i64 or cast to __int64. + // result of 32-bit shift implicitly converted to 64 bits (was 64-bit shift intended?) + uint_fast64_t mask = 1; + mask = mask << (index % 64); + if (value) { + bucket_array[bucket] |= mask; + } else { + bucket_array[bucket] &= ~mask; + } + } + + bool get(const uint_fast64_t index) { + uint_fast64_t bucket = index / 64; + // Taking the step with mask is crucial as we NEED a 64bit shift, not a 32bit one. + // MSVC: C4334, use 1i64 or cast to __int64. + // result of 32-bit shift implicitly converted to 64 bits (was 64-bit shift intended?) + uint_fast64_t mask = 1; + mask = mask << (index % 64); + return ((bucket_array[bucket] & mask) == mask); + } + + // Operators + BitVector operator &(BitVector const &bv) { + uint_fast64_t minSize = (bv.bucket_count < this->bucket_count) ? bv.bucket_count : this->bucket_count; + BitVector result(minSize * 64); + for (uint_fast64_t i = 0; i < minSize; ++i) { + result.bucket_array[i] = this->bucket_array[i] & bv.bucket_array[i]; + } + + return result; + } + BitVector operator |(BitVector const &bv) { + uint_fast64_t minSize = (bv.bucket_count < this->bucket_count) ? bv.bucket_count : this->bucket_count; + BitVector result(minSize * 64); + for (uint_fast64_t i = 0; i < minSize; ++i) { + result.bucket_array[i] = this->bucket_array[i] | bv.bucket_array[i]; + } + + return result; + } + + BitVector operator ^(BitVector const &bv) { + uint_fast64_t minSize = (bv.bucket_count < this->bucket_count) ? bv.bucket_count : this->bucket_count; + BitVector result(minSize * 64); + for (uint_fast64_t i = 0; i < minSize; ++i) { + result.bucket_array[i] = this->bucket_array[i] ^ bv.bucket_array[i]; + } + + return result; + } + + BitVector operator ~() { + BitVector result(this->bucket_count * 64); + for (uint_fast64_t i = 0; i < this->bucket_count; ++i) { + result.bucket_array[i] = ~this->bucket_array[i]; + } + + return result; + } + private: + uint_fast64_t bucket_count; + + /*! Array containing the boolean bits for each node, 64bits/64nodes per element */ + uint_fast64_t* bucket_array; + +}; + +} // namespace vector + +} // namespace mrmc + +#endif // MRMC_SPARSE_STATIC_SPARSE_MATRIX_H_ diff --git a/test/eigen/sparse_matrix_test.cpp b/test/eigen/sparse_matrix_test.cpp new file mode 100644 index 000000000..25928d06a --- /dev/null +++ b/test/eigen/sparse_matrix_test.cpp @@ -0,0 +1,79 @@ +#include "gtest/gtest.h" + +#include "Eigen/Sparse" +#include "src/exceptions/invalid_argument.h" +#include "boost/integer/integer_mask.hpp" + +TEST(EigenSparseMatrixTest, BasicReadWriteTest) { + // 25 rows, 50 non zero entries + Eigen::SparseMatrix esm(25, 25); + + int values[50] = { + 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, + 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, + 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, + 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, + 40, 41, 42, 43, 44, 45, 46, 47, 48, 49 + }; + 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 */ + 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, + 24, 24, 24, 24, 24, 24, 24, 24, 24, 24 /* second to last row */ + }; + int position_col[50] = { + 1, 3, 4, 5, 6, 7, 8, 9, 10, + 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, + 21, 22, 23, 24, 25, /* first row empty, one full row a 25 */ + 1, 25, /* one empty row, then first and last column */ + 16, 17, 18, 19, /* a few empty rows, middle columns */ + 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, + 14, 15, 16, 17, 18, 19, 20, 21, 22, 23 /* second to last row */ + }; + + typedef Eigen::Triplet ETd; + std::vector tripletList; + tripletList.reserve(50); + + for (int i = 0; i < 50; ++i) { + ASSERT_NO_THROW(tripletList.push_back(ETd(position_row[i] - 1, position_col[i] - 1, values[i]))); + } + + ASSERT_NO_THROW(esm.setFromTriplets(tripletList.begin(), tripletList.end())); + + for (int i = 0; i < 50; ++i) { + Eigen::SparseMatrix::InnerIterator it(esm, position_row[i] - 1); + ASSERT_EQ(values[i], esm.coeff(position_row[i] - 1, position_col[i] - 1)); + } + + // test for a few of the empty rows + for (int row = 15; row < 24; ++row) { + for (int col = 1; col <= 25; ++col) { + ASSERT_EQ(0, esm.coeff(row - 1, col - 1)); + } + } +} + +TEST(EigenSparseMatrixTest, BoundaryTest) { + Eigen::SparseMatrix esm(10, 10); + esm.reserve(100); + int values[100]; + for (uint_fast32_t i = 0; i < 100; ++i) { + values[i] = i; + } + + for (uint_fast32_t row = 0; row < 10; ++row) { + for (uint_fast32_t col = 0; col < 10; ++col) { + ASSERT_NO_THROW(esm.insert(row, col) = values[row * 10 + col]); + } + } + + for (uint_fast32_t row = 0; row < 10; ++row) { + for (uint_fast32_t col = 0; col < 10; ++col) { + ASSERT_EQ(values[row * 10 + col], esm.coeff(row, col)); + } + } +} \ No newline at end of file diff --git a/test/vector/bitvector_test.cpp b/test/vector/bitvector_test.cpp new file mode 100644 index 000000000..d25183b02 --- /dev/null +++ b/test/vector/bitvector_test.cpp @@ -0,0 +1,114 @@ +#include "gtest/gtest.h" +#include "src/vector/bitvector.h" +#include "src/exceptions/invalid_argument.h" + +TEST(BitVectorTest, GetSetTest) { + mrmc::vector::BitVector *bv = NULL; + ASSERT_NO_THROW(bv = new mrmc::vector::BitVector(32)); + + for (int i = 0; i < 32; ++i) { + bv->set(i, i % 2 == 0); + } + + for (int i = 0; i < 32; ++i) { + ASSERT_EQ(bv->get(i), i % 2 == 0); + } + delete bv; +} + +TEST(BitVectorTest, InitialZeroTest) { + mrmc::vector::BitVector bvA(32); + + for (int i = 0; i < 32; ++i) { + ASSERT_FALSE(bvA.get(i)); + } +} + +TEST(BitVectorTest, ResizeTest) { + mrmc::vector::BitVector bvA(32); + + for (int i = 0; i < 32; ++i) { + bvA.set(i, true); + } + + bvA.resize(70); + + for (int i = 0; i < 32; ++i) { + ASSERT_TRUE(bvA.get(i)); + } + + bool result; + for (int i = 32; i < 70; ++i) { + result = true; + ASSERT_NO_THROW(result = bvA.get(i)); + ASSERT_FALSE(result); + } +} + +TEST(BitVectorTest, OperatorNotTest) { + mrmc::vector::BitVector bvA(32); + mrmc::vector::BitVector bvB(32); + + for (int i = 0; i < 32; ++i) { + bvA.set(i, i % 2 == 0); + bvB.set(i, i % 2 == 1); + } + + mrmc::vector::BitVector bvN = ~bvB; + + for (int i = 0; i < 32; ++i) { + ASSERT_EQ(bvA.get(i), bvN.get(i)); + } +} + +TEST(BitVectorTest, OperatorAndTest) { + mrmc::vector::BitVector bvA(32); + mrmc::vector::BitVector bvB(32); + + for (int i = 0; i < 32; ++i) { + bvA.set(i, i % 2 == 0); + bvB.set(i, i % 2 == 1); + } + + mrmc::vector::BitVector bvN = bvA & bvB; + + for (int i = 0; i < 32; ++i) { + ASSERT_FALSE(bvN.get(i)); + } +} + +TEST(BitVectorTest, OperatorOrTest) { + mrmc::vector::BitVector bvA(32); + mrmc::vector::BitVector bvB(32); + + for (int i = 0; i < 32; ++i) { + bvA.set(i, i % 2 == 0); + bvB.set(i, i % 2 == 1); + } + + mrmc::vector::BitVector bvN = bvA | bvB; + + for (int i = 0; i < 32; ++i) { + ASSERT_TRUE(bvN.get(i)); + } +} + +TEST(BitVectorTest, OperatorXorTest) { + mrmc::vector::BitVector bvA(32); + mrmc::vector::BitVector bvB(32); + + for (int i = 0; i < 32; ++i) { + bvA.set(i, true); + bvB.set(i, i % 2 == 1); + } + + mrmc::vector::BitVector bvN = bvA ^ bvB; + mrmc::vector::BitVector bvO = ~bvB; + mrmc::vector::BitVector bvP = bvA ^ bvA; + + for (int i = 0; i < 32; ++i) { + ASSERT_EQ(bvN.get(i), bvO.get(i)); + // A XOR A = 0 + ASSERT_FALSE(bvP.get(i)); + } +} \ No newline at end of file