Browse Source

Added BitVector, a simple boolean array

Edited AtomicProposition to make use of BitVector
Edited mrmc-cpp.cpp, now comes with speed tests for some tasks
Edited static_sparse_matrix.h, added a copy constructor, added template  conversion functions, applied template parameters to all relevant functions/variables
Removed eigen_sparse_additions.h from internal files
tempestpy_adaptions
PBerger 12 years ago
parent
commit
6d384967fc
  1. 41
      src/dtmc/atomic_proposition.h
  2. 95
      src/mrmc-cpp.cpp
  3. 264
      src/sparse/static_sparse_matrix.h
  4. 152
      src/vector/bitvector.h
  5. 79
      test/eigen/sparse_matrix_test.cpp
  6. 114
      test/vector/bitvector_test.cpp

41
src/dtmc/atomic_proposition.h

@ -8,10 +8,7 @@
#include <pantheios/pantheios.hpp>
#include <pantheios/inserters/integer.hpp>
#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<int>(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<int>(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

95
src/mrmc-cpp.cpp

@ -15,7 +15,8 @@
#include <iostream>
#include <cstdio>
#include <ctime>
/* PlatformSTL Header Files */
#include <platformstl/performance/performance_counter.hpp>
#include <pantheios/pantheios.hpp>
#include <pantheios/backends/bec.file.h>
@ -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<double> *ssm;
mrmc::sparse::StaticSparseMatrix<double> *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<int> *ssm = new mrmc::sparse::StaticSparseMatrix<int>(10, 10);
ssmBack = new mrmc::sparse::StaticSparseMatrix<double>(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;
}
free(iArray);
}
printf("\n");
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();
delete ssm;
*/
std::cout << "Cpp-style: " << (unsigned long long)counter.get_microseconds() << std::endl;
return 0;
}

264
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<T> &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<T> &).
\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<typename _Scalar, int _Options, typename _Index>
void initialize(const Eigen::SparseMatrix<_Scalar, _Options, _Index> &eigen_sparse_matrix) {
template<int _Options, typename _Index>
void initialize(const Eigen::SparseMatrix<T, _Options, _Index> &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<T> toEigenSparseMatrix() {
Eigen::SparseMatrix<T> mat(row_count, row_count);
Eigen::SparseMatrix<T, Eigen::RowMajor, int_fast32_t>* toEigenSparseMatrix() {
int_fast32_t eigenRows = static_cast<int_fast32_t>(row_count);
Eigen::SparseMatrix<T, Eigen::RowMajor, int_fast32_t>* mat = new Eigen::SparseMatrix<T, Eigen::RowMajor, int_fast32_t>(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<int> IntTriplet;
# ifdef MRMC_USE_TRIPLETCONVERT
typedef Eigen::Triplet<T> IntTriplet;
std::vector<IntTriplet> 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<int_fast32_t>(mat->outerSize()), static_cast<int_fast32_t>((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<T> 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<T>)
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<T> 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_fast64_t getNonZeroEntryCount() const {
return non_zero_entry_count;
}
uint_fast32_t eigen_non_zero_entries = mrmc::sparse::getEigenSparseMatrixCorrectNonZeroEntryCount(eigen_sparse_matrix);
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<T>::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<typename _Scalar>
_Scalar constGetZero() const {
return _Scalar(0);
}
template <>
int_fast32_t constGetZero() const {
return 0;
}
template <>
double constGetZero() const {
return 0.0;
}
template<typename _Scalar>
_Scalar constGetOne() const {
return _Scalar(1);
}
template<>
int_fast32_t constGetOne() const {
return 1;
}
template<>
double constGetOne() const {
return 1.0;
}
template <typename _Scalar, typename _Index>
bool isEigenRowMajor(Eigen::SparseMatrix<_Scalar, Eigen::RowMajor, _Index>) {
return true;
@ -477,21 +533,21 @@ class StaticSparseMatrix {
}
template<typename _Scalar, int _Options, typename _Index>
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<uint_fast32_t>(entryCount - diag_non_zeros);
return static_cast<_Index>(entryCount - diag_non_zeros);
}
};

152
src/vector/bitvector.h

@ -0,0 +1,152 @@
#ifndef MRMC_VECTOR_BITVECTOR_H_
#define MRMC_VECTOR_BITVECTOR_H_
#include <exception>
#include <new>
#include <cmath>
#include "boost/integer/integer_mask.hpp"
#include <pantheios/pantheios.hpp>
#include <pantheios/inserters/integer.hpp>
#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_

79
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<int, 1> 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<double> ETd;
std::vector<ETd> 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<int, 1>::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<int, 1> 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));
}
}
}

114
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));
}
}
Loading…
Cancel
Save