From c30d6d307e5060ef8502b22f7fbab9d231ce4118 Mon Sep 17 00:00:00 2001 From: masawei Date: Mon, 18 Nov 2013 23:27:45 +0100 Subject: [PATCH] Figured out how to explicitly instantiate templates. But got bitten by std::vector as it is specialized and uses bitsets (i.e. integers) internally. Less memory but at the cost of 'oh, sorry std::vector does not return a bool&'. That again seems to be a problem for the SparseMatrix instatiation since for instance getValue returns a T&. On the one hand I don't quite know why this was never an issue before and on the other hand it prevents successful compilation. So there are different ways to settle this: - Specialize SparseMatix for bool -> possibly lots of code, but might be the best solution - Write a wrapper for std::vector that uses chars instead of booleans - Dont't use SparseMatrix Next up: Figure out the best solution for this and implement it. Former-commit-id: 83b9cfd06ead5399502f3836bfe243fcfacacab1 --- src/storage/SparseMatrix.cpp | 1040 +++++++++++++++++ src/storage/SparseMatrix.h | 950 ++------------- .../StronglyConnectedComponentDecomposition.h | 3 +- src/storm.cpp | 2 +- 4 files changed, 1132 insertions(+), 863 deletions(-) create mode 100644 src/storage/SparseMatrix.cpp diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp new file mode 100644 index 000000000..6cc37eca9 --- /dev/null +++ b/src/storage/SparseMatrix.cpp @@ -0,0 +1,1040 @@ +/* + * SparseMatrix.cpp + * + * Created on: Nov 18, 2013 + * Author: Manuel Sascha Weiand + */ + +#include "src/storage/SparseMatrix.h" + +namespace storm { +namespace storage { + + // Functions of the nested ConstIterator class. + + template + SparseMatrix::ConstIterator::ConstIterator(T const* valuePtr, uint_fast64_t const* columnPtr) : valuePtr(valuePtr), columnPtr(columnPtr) { + // Intentionally left empty. + } + + template + typename SparseMatrix::ConstIterator& SparseMatrix::ConstIterator::operator++() { + this->valuePtr++; + this->columnPtr++; + return *this; + } + + template + typename SparseMatrix::ConstIterator& SparseMatrix::ConstIterator::operator*() { + return *this; + } + + template + bool SparseMatrix::ConstIterator::operator!=(ConstIterator const& other) const { + return this->valuePtr != other.valuePtr; + } + + template + typename SparseMatrix::ConstIterator& SparseMatrix::ConstIterator::operator=(ConstIterator const& other) { + this->valuePtr = other.valuePtr; + this->columnPtr = other.columnPtr; + return *this; + } + + template + uint_fast64_t SparseMatrix::ConstIterator::column() const { + return *columnPtr; + } + + template + T const& SparseMatrix::ConstIterator::value() const { + return *valuePtr; + } + + + // Functions of the nested Rows class. + + template + SparseMatrix::Rows::Rows(T const* valuePtr, uint_fast64_t const* columnPtr, uint_fast64_t entryCount) : valuePtr(valuePtr), columnPtr(columnPtr), entryCount(entryCount) { + // Intentionally left empty. + } + + template + typename SparseMatrix::ConstIterator SparseMatrix::Rows::begin() const { + return ConstIterator(valuePtr, columnPtr); + } + + template + typename SparseMatrix::ConstIterator SparseMatrix::Rows::end() const { + return ConstIterator(valuePtr + entryCount, columnPtr + entryCount); + } + + + // Functions of the nested ConstRowsIterator class. + + template + SparseMatrix::ConstRowIterator::ConstRowIterator(T const* startValuePtr, uint_fast64_t const* startColumnPtr, uint_fast64_t const* rowPtr) : startValuePtr(startValuePtr), startColumnPtr(startColumnPtr), rowPtr(rowPtr) { + // Intentionally left empty. + } + + template + typename SparseMatrix::ConstRowIterator& SparseMatrix::ConstRowIterator::operator++() { + ++rowPtr; + return *this; + } + + template + bool SparseMatrix::ConstRowIterator::operator!=(ConstRowIterator const& other) const { + return this->rowPtr != other.rowPtr; + } + + template + typename SparseMatrix::ConstIterator SparseMatrix::ConstRowIterator::begin() const { + return ConstIterator(startValuePtr + *rowPtr, startColumnPtr + *rowPtr); + } + + template + typename SparseMatrix::ConstIterator SparseMatrix::ConstRowIterator::end() const { + return ConstIterator(startValuePtr + *(rowPtr + 1), startColumnPtr + *(rowPtr + 1)); + } + + + // Functions of the SparseMatrix class. + + template + SparseMatrix::SparseMatrix(uint_fast64_t rows, uint_fast64_t cols) : rowCount(rows), colCount(cols), + nonZeroEntryCount(0), internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { + // Intentionally left empty. + } + + template + SparseMatrix::SparseMatrix(uint_fast64_t size) + : rowCount(size), colCount(size), nonZeroEntryCount(0), + internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { + // Intentionally left empty. + } + + template + SparseMatrix::SparseMatrix(SparseMatrix&& other) + : rowCount(other.rowCount), colCount(other.colCount), nonZeroEntryCount(other.nonZeroEntryCount), + valueStorage(std::move(other.valueStorage)), columnIndications(std::move(other.columnIndications)), + rowIndications(std::move(other.rowIndications)), internalStatus(other.internalStatus), + currentSize(other.currentSize), lastRow(other.lastRow) { + // Now update the source matrix + other.rowCount = 0; + other.colCount = 0; + other.nonZeroEntryCount = 0; + other.internalStatus = MatrixStatus::Error; + other.currentSize = 0; + other.lastRow = 0; + } + + template + SparseMatrix::SparseMatrix(const SparseMatrix & other) + : rowCount(other.rowCount), colCount(other.colCount), nonZeroEntryCount(other.nonZeroEntryCount), + valueStorage(other.valueStorage), columnIndications(other.columnIndications), + rowIndications(other.rowIndications), internalStatus(other.internalStatus), + currentSize(other.currentSize), lastRow(other.lastRow) { + } + + template + SparseMatrix::SparseMatrix(uint_fast64_t rowCount, uint_fast64_t colCount, uint_fast64_t nonZeroEntryCount, + std::vector&& rowIndications, + std::vector&& columnIndications, std::vector&& values) + : rowCount(rowCount), colCount(colCount), nonZeroEntryCount(nonZeroEntryCount), + valueStorage(values), columnIndications(columnIndications), + rowIndications(rowIndications), internalStatus(MatrixStatus::Initialized), + currentSize(0), lastRow(0) { + // Intentionally left empty. + } + + template + storm::storage::SparseMatrix& SparseMatrix::operator=(SparseMatrix const& other) { + this->rowCount = other.rowCount; + this->colCount = other.colCount; + this->nonZeroEntryCount = other.nonZeroEntryCount; + + this->valueStorage = other.valueStorage; + this->columnIndications = other.columnIndications; + this->rowIndications = other.rowIndications; + + this->internalStatus = other.internalStatus; + this->currentSize = other.currentSize; + this->lastRow = other.lastRow; + + return *this; + } + + template + void SparseMatrix::initialize(uint_fast64_t nonZeroEntries) { + // Check whether initializing the matrix is safe. + if (internalStatus != MatrixStatus::UnInitialized) { + triggerErrorState(); + LOG4CPLUS_ERROR(logger, "Trying to initialize matrix that is not uninitialized."); + throw storm::exceptions::InvalidStateException("Trying to initialize matrix that is not uninitialized."); + } else if ((rowCount * colCount) < nonZeroEntries) { + triggerErrorState(); + LOG4CPLUS_ERROR(logger, "Trying to initialize a matrix with more non-zero entries than there can be."); + throw storm::exceptions::InvalidArgumentException("Trying to initialize a matrix with more non-zero entries than there can be."); + } else { + // If it is safe, initialize necessary members and prepare the + // internal storage. + nonZeroEntryCount = nonZeroEntries; + lastRow = 0; + + if (!prepareInternalStorage()) { + triggerErrorState(); + throw std::bad_alloc(); + } else { + setState(MatrixStatus::Initialized); + } + } + } + + template + void SparseMatrix::addNextValue(const uint_fast64_t row, const uint_fast64_t col, T const& value) { + // Check whether the given row and column positions are valid and throw + // error otherwise. + if ((row > rowCount) || (col > colCount)) { + triggerErrorState(); + LOG4CPLUS_ERROR(logger, "Trying to add a value at illegal position (" << row << ", " << col << ") in matrix of size (" << rowCount << ", " << colCount << ")."); + throw storm::exceptions::OutOfRangeException() << "Trying to add a value at illegal position (" << row << ", " << col << ") in matrix of size (" << rowCount << ", " << colCount << ")."; + } + + // 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; + } + lastRow = row; + } + + // Finally, set the element and increase the current size. + valueStorage[currentSize] = value; + columnIndications[currentSize] = col; + + ++currentSize; + } + + template + void SparseMatrix::insertNextValue(const uint_fast64_t row, const uint_fast64_t col, T const& value, bool pushRowIndication) { + // Check whether the given row and column positions are valid and throw + // error otherwise. + if (row < lastRow) { + triggerErrorState(); + LOG4CPLUS_ERROR(logger, "Trying to insert a value at illegal position (" << row << ", " << col << ")."); + throw storm::exceptions::OutOfRangeException() << "Trying to insert a value at illegal position (" << row << ", " << col << ")."; + } + + // If we switched to another row, we have to adjust the missing entries in the rowIndications array. + if (row != lastRow) { + for (uint_fast64_t i = lastRow + 1; i <= row; ++i) { + if (pushRowIndication) { + rowIndications.push_back(currentSize); + } else { + rowIndications[i] = currentSize; + } + } + rowCount = row + 1; + lastRow = row; + } + + // Finally, set the element and increase the current size. + valueStorage.push_back(value); + columnIndications.push_back(col); + ++nonZeroEntryCount; + ++currentSize; + + // Check that we also have the correct number of columns. + colCount = std::max(colCount, col + 1); + } + + template + void SparseMatrix::insertEmptyRow(bool pushRowIndication) { + if (pushRowIndication) { + rowIndications.push_back(currentSize); + } else { + rowIndications[lastRow + 1] = currentSize; + } + + ++rowCount; + ++lastRow; + } + + template + void SparseMatrix::finalize(bool pushSentinelElement) { + // Check whether it's safe to finalize the matrix and throw error otherwise. + if (!isInitialized()) { + triggerErrorState(); + LOG4CPLUS_ERROR(logger, "Trying to finalize an uninitialized matrix."); + throw storm::exceptions::InvalidStateException("Trying to finalize an uninitialized matrix."); + } else if (currentSize != nonZeroEntryCount) { + triggerErrorState(); + LOG4CPLUS_ERROR(logger, "Trying to finalize a matrix that was initialized with more non-zero entries than given (expected " << nonZeroEntryCount << " but got " << currentSize << " instead)"); + throw storm::exceptions::InvalidStateException() << "Trying to finalize a matrix that was initialized with more non-zero entries than given (expected " << nonZeroEntryCount << " but got " << currentSize << " instead)."; + } else { + // Fill in the missing entries in the row_indications array. + // (Can happen because of empty rows at the end.) + for (uint_fast64_t i = lastRow + 1; i < rowCount; ++i) { + rowIndications[i] = currentSize; + } + + // 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 rowIndications[i] and rowIndications[i + 1], also for the + // first and last row. + if (pushSentinelElement) { + rowIndications.push_back(nonZeroEntryCount); + } else { + rowIndications[rowCount] = nonZeroEntryCount; + } + + setState(MatrixStatus::ReadReady); + } + } + + template + inline bool SparseMatrix::getValue(uint_fast64_t row, uint_fast64_t col, T* const target) const { + // Check for illegal access indices. + if ((row > rowCount) || (col > colCount)) { + LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ")."); + throw storm::exceptions::OutOfRangeException("Trying to read a value from illegal position."); + return false; + } + + // 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 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; + } + ++rowStart; + } + + // Set 0 as the content and return false in case the element was not found. + *target = 0; + return false; + } + + template + inline T& SparseMatrix::getValue(uint_fast64_t row, uint_fast64_t col) { + // Check for illegal access indices. + if ((row > rowCount) || (col > colCount)) { + LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ")."); + throw storm::exceptions::OutOfRangeException("Trying to read a value from illegal position."); + } + + // 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, return it. + if (columnIndications[rowStart] == col) { + return valueStorage[rowStart]; + } + + // 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; + } + ++rowStart; + } + + throw storm::exceptions::InvalidArgumentException("Trying to get a reference to a non-existant value."); + } + + template + uint_fast64_t SparseMatrix::getRowCount() const { + return rowCount; + } + + template + uint_fast64_t SparseMatrix::getColumnCount() const { + return colCount; + } + + template + bool SparseMatrix::isReadReady() { + return (internalStatus == MatrixStatus::ReadReady); + } + + template + bool SparseMatrix::isInitialized() { + return (internalStatus == MatrixStatus::Initialized || internalStatus == MatrixStatus::ReadReady); + } + + template + typename SparseMatrix::MatrixStatus SparseMatrix::getState() { + return internalStatus; + } + + template + bool SparseMatrix::hasError() const { + return (internalStatus == MatrixStatus::Error); + } + + template + uint_fast64_t SparseMatrix::getNonZeroEntryCount() const { + return nonZeroEntryCount; + } + + template + bool SparseMatrix::makeRowsAbsorbing(storm::storage::BitVector const& rows) { + bool result = true; + for (auto row : rows) { + result &= makeRowAbsorbing(row, row); + } + return result; + } + + template + bool SparseMatrix::makeRowsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices) { + bool result = true; + for (auto rowGroup : rowGroupConstraint) { + for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) { + result &= makeRowAbsorbing(row, rowGroup); + } + } + return result; + } + + template + bool SparseMatrix::makeRowAbsorbing(const uint_fast64_t row, const uint_fast64_t column) { + // Check whether the accessed state exists. + if (row > rowCount) { + LOG4CPLUS_ERROR(logger, "Trying to make an illegal row " << row << " absorbing."); + throw storm::exceptions::OutOfRangeException() << "Trying to make an illegal row " << row << " absorbing."; + return false; + } + + // Iterate over the elements in the row that are not on the diagonal + // and set them to zero. + uint_fast64_t rowStart = rowIndications[row]; + uint_fast64_t rowEnd = rowIndications[row + 1]; + + // If the row has no elements in it, we cannot make it absorbing, because we would need to + // move all elements in the vector of nonzeros otherwise. + if (rowStart >= rowEnd) { + LOG4CPLUS_ERROR(logger, "Cannot make row " << row << " absorbing, because there is no entry in this row."); + throw storm::exceptions::InvalidStateException() << "Cannot make row " << row << " absorbing, because there is no entry in this row."; + } + + // If there is at least one nonzero entry in this row, we can just set it to one, modify its + // column indication to the one given by the parameter and set all subsequent elements of this + // row to zero. + valueStorage[rowStart] = storm::utility::constGetOne(); + columnIndications[rowStart] = column; + for (uint_fast64_t index = rowStart + 1; index < rowEnd; ++index) { + valueStorage[index] = storm::utility::constGetZero(); + columnIndications[index] = 0; + } + + return true; + } + + template + T SparseMatrix::getConstrainedRowSum(uint_fast64_t row, storm::storage::BitVector const& constraint) const { + T result(0); + for (uint_fast64_t i = rowIndications[row]; i < rowIndications[row + 1]; ++i) { + if (constraint.get(columnIndications[i])) { + result += valueStorage[i]; + } + } + return result; + } + + template + std::vector SparseMatrix::getConstrainedRowSumVector(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint) const { + std::vector result(rowConstraint.getNumberOfSetBits()); + uint_fast64_t currentRowCount = 0; + for (auto row : rowConstraint) { + result[currentRowCount++] = getConstrainedRowSum(row, columnConstraint); + } + return result; + } + + template + std::vector SparseMatrix::getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices, storm::storage::BitVector const& columnConstraint, uint_fast64_t numberOfRows) const { + std::vector result(numberOfRows); + uint_fast64_t currentRowCount = 0; + for (auto rowGroup : rowGroupConstraint) { + for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) { + result[currentRowCount++] = getConstrainedRowSum(row, columnConstraint); + } + } + return result; + } + + template + SparseMatrix SparseMatrix::getSubmatrix(storm::storage::BitVector const& constraint) const { + LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix with " << constraint.getNumberOfSetBits() << " rows."); + + // Check for valid constraint. + if (constraint.getNumberOfSetBits() == 0) { + LOG4CPLUS_ERROR(logger, "Trying to create a sub-matrix of size 0."); + throw storm::exceptions::InvalidArgumentException("Trying to create a sub-matrix of size 0."); + } + + // First, we need to determine the number of non-zero entries of the + // sub-matrix. + uint_fast64_t subNonZeroEntries = 0; + for (auto rowIndex : constraint) { + for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) { + if (constraint.get(columnIndications[i])) { + ++subNonZeroEntries; + } + } + } + + // Create and initialize resulting matrix. + SparseMatrix result(constraint.getNumberOfSetBits()); + result.initialize(subNonZeroEntries); + + // Create a temporary vecotr that stores for each index whose bit is set + // to true the number of bits that were set before that particular index. + std::vector bitsSetBeforeIndex; + bitsSetBeforeIndex.reserve(colCount); + + // Compute the information to fill this vector. + uint_fast64_t lastIndex = 0; + uint_fast64_t currentNumberOfSetBits = 0; + for (auto index : constraint) { + while (lastIndex <= index) { + bitsSetBeforeIndex.push_back(currentNumberOfSetBits); + ++lastIndex; + } + ++currentNumberOfSetBits; + } + + // Copy over selected entries. + uint_fast64_t rowCount = 0; + for (auto rowIndex : constraint) { + for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) { + if (constraint.get(columnIndications[i])) { + result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[i]], valueStorage[i]); + } + } + + ++rowCount; + } + + // Finalize sub-matrix and return result. + result.finalize(); + LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); + return result; + } + + template + SparseMatrix SparseMatrix::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices) const { + LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size)."); + + // First, we need to determine the number of non-zero entries and the number of rows of the sub-matrix. + uint_fast64_t subNonZeroEntries = 0; + uint_fast64_t subRowCount = 0; + for (auto index : rowGroupConstraint) { + subRowCount += rowGroupIndices[index + 1] - rowGroupIndices[index]; + for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { + for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { + if (rowGroupConstraint.get(columnIndications[j])) { + ++subNonZeroEntries; + } + } + } + } + + LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << rowGroupConstraint.getNumberOfSetBits() << "."); + + // Create and initialize resulting matrix. + SparseMatrix result(subRowCount, rowGroupConstraint.getNumberOfSetBits()); + result.initialize(subNonZeroEntries); + + // Create a temporary vector that stores for each index whose bit is set + // to true the number of bits that were set before that particular index. + std::vector bitsSetBeforeIndex; + bitsSetBeforeIndex.reserve(colCount); + + // Compute the information to fill this vector. + uint_fast64_t lastIndex = 0; + uint_fast64_t currentNumberOfSetBits = 0; + for (auto index : rowGroupConstraint) { + while (lastIndex <= index) { + bitsSetBeforeIndex.push_back(currentNumberOfSetBits); + ++lastIndex; + } + ++currentNumberOfSetBits; + } + + // Copy over selected entries. + uint_fast64_t rowCount = 0; + for (auto index : rowGroupConstraint) { + for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { + for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { + if (rowGroupConstraint.get(columnIndications[j])) { + result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]); + } + } + ++rowCount; + } + } + + // Finalize sub-matrix and return result. + result.finalize(); + LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); + return result; + } + + template + SparseMatrix SparseMatrix::getSubmatrix(std::vector const& rowGroupToRowIndexMapping, std::vector const& rowGroupIndices, bool insertDiagonalEntries) const { + LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size)."); + + // First, we need to count how many non-zero entries the resulting matrix will have and reserve space for diagonal + // entries. + uint_fast64_t subNonZeroEntries = 0; + for (uint_fast64_t rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) { + // Determine which row we need to select from the current row group. + uint_fast64_t rowToCopy = rowGroupIndices[rowGroupIndex] + rowGroupToRowIndexMapping[rowGroupIndex]; + + // Iterate through that row and count the number of slots we have to reserve for copying. + bool foundDiagonalElement = false; + for (uint_fast64_t i = rowIndications[rowToCopy], rowEnd = rowIndications[rowToCopy + 1]; i < rowEnd; ++i) { + if (columnIndications[i] == rowGroupIndex) { + foundDiagonalElement = true; + } + ++subNonZeroEntries; + } + if (insertDiagonalEntries && !foundDiagonalElement) { + ++subNonZeroEntries; + } + } + + LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << (rowGroupIndices.size() - 1) << "x" << colCount << " with " << subNonZeroEntries << " non-zero elements."); + + // Now create the matrix to be returned with the appropriate size. + SparseMatrix submatrix(rowGroupIndices.size() - 1, colCount); + submatrix.initialize(subNonZeroEntries); + + // Copy over the selected lines from the source matrix. + for (uint_fast64_t rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) { + // Determine which row we need to select from the current row group. + uint_fast64_t rowToCopy = rowGroupIndices[rowGroupIndex] + rowGroupToRowIndexMapping[rowGroupIndex]; + + // Iterate through that row and copy the entries. This also inserts a zero element on the diagonal if there + // is no entry yet. + bool insertedDiagonalElement = false; + for (uint_fast64_t i = rowIndications[rowToCopy], rowEnd = rowIndications[rowToCopy + 1]; i < rowEnd; ++i) { + if (columnIndications[i] == rowGroupIndex) { + insertedDiagonalElement = true; + } else if (insertDiagonalEntries && !insertedDiagonalElement && columnIndications[i] > rowGroupIndex) { + submatrix.addNextValue(rowGroupIndex, rowGroupIndex, storm::utility::constGetZero()); + insertedDiagonalElement = true; + } + submatrix.addNextValue(rowGroupIndex, columnIndications[i], valueStorage[i]); + } + if (insertDiagonalEntries && !insertedDiagonalElement) { + submatrix.addNextValue(rowGroupIndex, rowGroupIndex, storm::utility::constGetZero()); + } + } + + // Finalize created matrix and return result. + submatrix.finalize(); + LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); + return submatrix; + } + + template + void SparseMatrix::convertToEquationSystem() { + invertDiagonal(); + negateAllNonDiagonalElements(); + } + + template + void SparseMatrix::invertDiagonal() { + // Check if the matrix is square, because only then it makes sense to perform this + // transformation. + if (this->getRowCount() != this->getColumnCount()) { + throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to be square!"; + } + + // Now iterate over all rows and set the diagonal elements to the inverted value. + // If there is a row without the diagonal element, an exception is thrown. + T one = storm::utility::constGetOne(); + bool foundDiagonalElement = false; + for (uint_fast64_t row = 0; row < rowCount; ++row) { + uint_fast64_t rowStart = rowIndications[row]; + uint_fast64_t rowEnd = rowIndications[row + 1]; + foundDiagonalElement = false; + while (rowStart < rowEnd) { + if (columnIndications[rowStart] == row) { + valueStorage[rowStart] = one - valueStorage[rowStart]; + foundDiagonalElement = true; + break; + } + ++rowStart; + } + + // Throw an exception if a row did not have an element on the diagonal. + if (!foundDiagonalElement) { + throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to contain all diagonal entries!"; + } + } + } + + template + void SparseMatrix::negateAllNonDiagonalElements() { + // Check if the matrix is square, because only then it makes sense to perform this + // transformation. + if (this->getRowCount() != this->getColumnCount()) { + throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to be square!"; + } + + // Iterate over all rows and negate all the elements that are not on the diagonal. + for (uint_fast64_t row = 0; row < rowCount; ++row) { + uint_fast64_t rowStart = rowIndications[row]; + uint_fast64_t rowEnd = rowIndications[row + 1]; + while (rowStart < rowEnd) { + if (columnIndications[rowStart] != row) { + valueStorage[rowStart] = - valueStorage[rowStart]; + } + ++rowStart; + } + } + } + + template + typename SparseMatrix::SparseJacobiDecomposition_t SparseMatrix::getJacobiDecomposition() const { + uint_fast64_t rowCount = this->getRowCount(); + uint_fast64_t colCount = this->getColumnCount(); + if (rowCount != colCount) { + throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::getJacobiDecomposition requires the Matrix to be square."; + } + storm::storage::SparseMatrix resultLU(*this); + storm::storage::SparseMatrix resultDinv(rowCount, colCount); + // no entries apart from those on the diagonal (rowCount) + resultDinv.initialize(rowCount); + + // constant 1 for diagonal inversion + T constOne = storm::utility::constGetOne(); + + // copy diagonal entries to other matrix + for (uint_fast64_t i = 0; i < rowCount; ++i) { + resultDinv.addNextValue(i, i, constOne / resultLU.getValue(i, i)); + resultLU.getValue(i, i) = storm::utility::constGetZero(); + } + resultDinv.finalize(); + + return std::make_pair(std::move(resultLU), std::move(resultDinv)); + } + + template + std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const { + // Prepare result. + std::vector result(rowCount, storm::utility::constGetZero()); + + // Iterate over all elements of the current matrix and either continue with the next element + // in case the given matrix does not have a non-zero element at this column position, or + // multiply the two entries and add the result to the corresponding position in the vector. + for (uint_fast64_t row = 0; row < rowCount && row < otherMatrix.rowCount; ++row) { + for (uint_fast64_t element = rowIndications[row], nextOtherElement = otherMatrix.rowIndications[row]; element < rowIndications[row + 1] && nextOtherElement < otherMatrix.rowIndications[row + 1]; ++element) { + if (columnIndications[element] < otherMatrix.columnIndications[nextOtherElement]) { + continue; + } else { + // If the precondition of this method (i.e. that the given matrix is a submatrix + // of the current one) was fulfilled, we know now that the two elements are in + // the same column, so we can multiply and add them to the row sum vector. + result[row] += otherMatrix.valueStorage[nextOtherElement] * valueStorage[element]; + ++nextOtherElement; + } + } + } + + return result; + } + + + template + void SparseMatrix::multiplyWithVector(std::vector const& vector, std::vector& result) const { +#ifdef STORM_HAVE_INTELTBB + tbb::parallel_for(tbb::blocked_range(0, result.size()), tbbHelper_MatrixRowVectorScalarProduct, std::vector, T>(this, &vector, &result)); +#else + ConstRowIterator rowIt = this->begin(); + + for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++rowIt) { + *resultIt = storm::utility::constGetZero(); + + for (auto elementIt = rowIt.begin(), elementIte = rowIt.end(); elementIt != elementIte; ++elementIt) { + *resultIt += elementIt.value() * vector[elementIt.column()]; + } + } +#endif + } + + template + uint_fast64_t SparseMatrix::getSizeInMemory() const { + uint_fast64_t size = sizeof(*this); + // Add value_storage size. + size += sizeof(T) * valueStorage.capacity(); + // Add column_indications size. + size += sizeof(uint_fast64_t) * columnIndications.capacity(); + // Add row_indications size. + size += sizeof(uint_fast64_t) * rowIndications.capacity(); + return size; + } + + template + typename SparseMatrix::Rows SparseMatrix::getRows(uint_fast64_t startRow, uint_fast64_t endRow) const { + return Rows(this->valueStorage.data() + this->rowIndications[startRow], this->columnIndications.data() + this->rowIndications[startRow], this->rowIndications[endRow + 1] - this->rowIndications[startRow]); + } + + template + typename SparseMatrix::Rows SparseMatrix::getRow(uint_fast64_t row) const { + return getRows(row, row); + } + + template + typename SparseMatrix::ConstRowIterator SparseMatrix::begin(uint_fast64_t initialRow) const { + return ConstRowIterator(this->valueStorage.data(), this->columnIndications.data(), this->rowIndications.data() + initialRow); + } + + template + typename SparseMatrix::ConstRowIterator SparseMatrix::end() const { + return ConstRowIterator(this->valueStorage.data(), this->columnIndications.data(), this->rowIndications.data() + rowCount); + } + + template + typename SparseMatrix::ConstRowIterator SparseMatrix::end(uint_fast64_t row) const { + return ConstRowIterator(this->valueStorage.data(), this->columnIndications.data(), this->rowIndications.data() + row + 1); + } + + template + typename SparseMatrix::ConstIndexIterator SparseMatrix::constColumnIteratorBegin(uint_fast64_t row) const { + return &(this->columnIndications[0]) + this->rowIndications[row]; + } + + template + typename SparseMatrix::ConstIndexIterator SparseMatrix::constColumnIteratorEnd() const { + return &(this->columnIndications[0]) + this->rowIndications[rowCount]; + } + + template + typename SparseMatrix::ConstIndexIterator SparseMatrix::constColumnIteratorEnd(uint_fast64_t row) const { + return &(this->columnIndications[0]) + this->rowIndications[row + 1]; + } + + template + typename SparseMatrix::ConstValueIterator SparseMatrix::constValueIteratorBegin(uint_fast64_t row) const { + return &(this->valueStorage[0]) + this->rowIndications[row]; + } + + template + typename SparseMatrix::ConstValueIterator SparseMatrix::constValueIteratorEnd() const { + return &(this->valueStorage[0]) + this->rowIndications[rowCount]; + } + + template + typename SparseMatrix::ConstValueIterator SparseMatrix::constValueIteratorEnd(uint_fast64_t row) const { + return &(this->valueStorage[0]) + this->rowIndications[row + 1]; + } + + template + T SparseMatrix::getRowSum(uint_fast64_t row) const { + T sum = storm::utility::constGetZero(); + for (auto it = this->constValueIteratorBegin(row), ite = this->constValueIteratorEnd(row); it != ite; ++it) { + sum += *it; + } + return sum; + } + + template + bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const { + // Check for matching sizes. + if (this->getRowCount() != matrix.getRowCount()) return false; + if (this->getColumnCount() != matrix.getColumnCount()) return false; + + // Check the subset property for all rows individually. + for (uint_fast64_t row = 0; row < this->getRowCount(); ++row) { + for (uint_fast64_t elem = rowIndications[row], elem2 = matrix.rowIndications[row]; elem < rowIndications[row + 1] && elem < matrix.rowIndications[row + 1]; ++elem) { + // Skip over all entries of the other matrix that are before the current entry in + // the current matrix. + while (elem2 < matrix.rowIndications[row + 1] && matrix.columnIndications[elem2] < columnIndications[elem]) { + ++elem2; + } + if (!(elem2 < matrix.rowIndications[row + 1]) || columnIndications[elem] != matrix.columnIndications[elem2]) { + return false; + } + } + } + return true; + } + + template + std::string SparseMatrix::toStringCompressed() const { + std::stringstream result; + result << rowIndications << std::endl; + result << columnIndications << std::endl; + result << valueStorage << std::endl; + return result.str(); + } + + template + std::string SparseMatrix::toString(std::vector const* rowGroupIndices) const { + std::stringstream result; + uint_fast64_t currentNondeterministicChoiceIndex = 0; + + // Print column numbers in header. + result << "\t\t"; + for (uint_fast64_t i = 0; i < colCount; ++i) { + result << i << "\t"; + } + result << std::endl; + + // Iterate over all rows. + for (uint_fast64_t i = 0; i < rowCount; ++i) { + uint_fast64_t nextIndex = rowIndications[i]; + + // If we need to group rows, print a dashed line in case we have moved to the next group of rows. + if (rowGroupIndices != nullptr) { + if (i == (*rowGroupIndices)[currentNondeterministicChoiceIndex]) { + if (i != 0) { + result << "\t(\t"; + for (uint_fast64_t j = 0; j < colCount - 2; ++j) { + result << "----"; + if (j == 1) { + result << "\t" << currentNondeterministicChoiceIndex << "\t"; + } + } + result << "\t)" << std::endl; + } + ++currentNondeterministicChoiceIndex; + } + } + + // Print the actual row. + result << i << "\t(\t"; + uint_fast64_t currentRealIndex = 0; + while (currentRealIndex < colCount) { + if (nextIndex < rowIndications[i + 1] && currentRealIndex == columnIndications[nextIndex]) { + result << valueStorage[nextIndex] << "\t"; + ++nextIndex; + } else { + result << "0\t"; + } + ++currentRealIndex; + } + result << "\t)\t" << i << std::endl; + } + + // Print column numbers in footer. + result << "\t\t"; + for (uint_fast64_t i = 0; i < colCount; ++i) { + result << i << "\t"; + } + result << std::endl; + + // Return final result. + return result.str(); + } + + template + std::size_t SparseMatrix::getHash() const { + std::size_t result = 0; + + boost::hash_combine(result, rowCount); + boost::hash_combine(result, colCount); + boost::hash_combine(result, nonZeroEntryCount); + boost::hash_combine(result, currentSize); + boost::hash_combine(result, lastRow); + boost::hash_combine(result, storm::utility::Hash::getHash(valueStorage)); + boost::hash_combine(result, storm::utility::Hash::getHash(columnIndications)); + boost::hash_combine(result, storm::utility::Hash::getHash(rowIndications)); + + return result; + } + + template + void SparseMatrix::triggerErrorState() { + setState(MatrixStatus::Error); + } + + template + void SparseMatrix::setState(const MatrixStatus new_state) { + internalStatus = (internalStatus == MatrixStatus::Error) ? internalStatus : new_state; + } + + template + bool SparseMatrix::prepareInternalStorage(bool initializeElements) { + if (initializeElements) { + // Set up the arrays for the elements that are not on the diagonal. + valueStorage.resize(nonZeroEntryCount, storm::utility::constGetZero()); + columnIndications.resize(nonZeroEntryCount, 0); + + // Set up the rowIndications vector and reserve one element more than there are rows in + // order to put a sentinel element at the end, which eases iteration process. + rowIndications.resize(rowCount + 1, 0); + + // Return whether all the allocations could be made without error. + return ((valueStorage.capacity() >= nonZeroEntryCount) && (columnIndications.capacity() >= nonZeroEntryCount) + && (rowIndications.capacity() >= (rowCount + 1))); + } else { + // If it was not requested to initialize the elements, we simply reserve the space. + valueStorage.reserve(nonZeroEntryCount); + columnIndications.reserve(nonZeroEntryCount); + rowIndications.reserve(rowCount + 1); + return true; + } + } + + // Explicit instantiations of specializations of this template here. + template class SparseMatrix; + template class SparseMatrix; + + // Functions of the tbbHelper_MatrixRowVectorScalarProduct friend class. + +#ifdef STORM_HAVE_INTELTBB + + template + tbbHelper_MatrixRowVectorScalarProduct::tbbHelper_MatrixRowVectorScalarProduct(M const* matrixA, V const* vectorX, V * resultVector) : matrixA(matrixA), vectorX(vectorX), resultVector(resultVector) {} + + template + void tbbHelper_MatrixRowVectorScalarProduct::operator() (const tbb::blocked_range& r) const { + for (uint_fast64_t row = r.begin(); row < r.end(); ++row) { + uint_fast64_t index = matrixA->rowIndications.at(row); + uint_fast64_t indexEnd = matrixA->rowIndications.at(row + 1); + + // Initialize the result to be 0. + T element = storm::utility::constGetZero(); + + for (; index != indexEnd; ++index) { + element += matrixA->valueStorage.at(index) * vectorX->at(matrixA->columnIndications.at(index)); + } + + // Write back to the result Vector + resultVector->at(row) = element; + } + } + + // Explicit instanciations of specializations of this template here. + template class tbbHelper_MatrixRowVectorScalarProduct, std::vector, double>; + +#endif + + +} // namespace storage +} // namespace storm + + + diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 6667fd82e..0bd275e73 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -99,20 +99,14 @@ public: * @param valuePtr A pointer to the value of the first element that is to be iterated over. * @param columnPtr A pointer to the column of the first element that is to be iterated over. */ - ConstIterator(T const* valuePtr, uint_fast64_t const* columnPtr) : valuePtr(valuePtr), columnPtr(columnPtr) { - // Intentionally left empty. - } + ConstIterator(T const* valuePtr, uint_fast64_t const* columnPtr); /*! * Moves the iterator to the next non-zero element. * * @return A reference to itself. */ - ConstIterator& operator++() { - ++valuePtr; - ++columnPtr; - return *this; - } + ConstIterator& operator++(); /*! * Dereferences the iterator by returning a reference to itself. This is needed for making use of the range-based @@ -120,29 +114,21 @@ public: * * @return A reference to itself. */ - ConstIterator& operator*() { - return *this; - } + ConstIterator& operator*(); /*! * Compares the two iterators for inequality. * * @return True iff the given iterator points to a different index as the current iterator. */ - bool operator!=(ConstIterator const& other) const { - return this->valuePtr != other.valuePtr; - } + bool operator!=(ConstIterator const& other) const; /*! * Assigns the position of the given iterator to the current iterator. * * @return A reference to itself. */ - ConstIterator& operator=(ConstIterator const& other) { - this->valuePtr = other.valuePtr; - this->columnPtr = other.columnPtr; - return *this; - } + ConstIterator& operator=(ConstIterator const& other); /*! * Retrieves the column that is associated with the current non-zero element to which this iterator @@ -150,18 +136,14 @@ public: * * @return The column of the current non-zero element to which this iterator points. */ - uint_fast64_t column() const { - return *columnPtr; - } + uint_fast64_t column() const; /*! * Retrieves the value of the current non-zero element to which this iterator points. * * @return The value of the current non-zero element to which this iterator points. */ - T const& value() const { - return *valuePtr; - } + T const& value() const; private: // A pointer to the value of the current non-zero element. @@ -183,27 +165,21 @@ public: * @param columnPtr A pointer to the column of the first non-zero element of the rows. * @param entryCount The number of non-zero elements of the rows. */ - Rows(T const* valuePtr, uint_fast64_t const* columnPtr, uint_fast64_t entryCount) : valuePtr(valuePtr), columnPtr(columnPtr), entryCount(entryCount) { - // Intentionally left empty. - } + Rows(T const* valuePtr, uint_fast64_t const* columnPtr, uint_fast64_t entryCount); /*! * Retrieves an iterator that points to the beginning of the rows. * * @return An iterator that points to the beginning of the rows. */ - ConstIterator begin() const { - return ConstIterator(valuePtr, columnPtr); - } + ConstIterator begin() const; /*! * Retrieves an iterator that points past the last element of the rows. * * @return An iterator that points past the last element of the rows. */ - ConstIterator end() const { - return ConstIterator(valuePtr + entryCount, columnPtr + entryCount); - } + ConstIterator end() const; private: // The pointer to the value of the first element. @@ -228,44 +204,33 @@ public: * @param startColumnPtr A pointer to the column of the first non-zero element of the matrix. * @param rowPtr A pointer to the index at which the first row begins. */ - ConstRowIterator(T const* startValuePtr, uint_fast64_t const* startColumnPtr, uint_fast64_t const* rowPtr) : startValuePtr(startValuePtr), startColumnPtr(startColumnPtr), rowPtr(rowPtr) { - // Intentionally left empty. - } + ConstRowIterator(T const* startValuePtr, uint_fast64_t const* startColumnPtr, uint_fast64_t const* rowPtr); /*! * This sets the iterator to point to the next row. * * @return A reference to itself. */ - ConstRowIterator& operator++() { - ++rowPtr; - return *this; - } + ConstRowIterator& operator++(); /*! * Compares the iterator to the given row iterator. */ - bool operator!=(ConstRowIterator const& other) const { - return this->rowPtr != other.rowPtr; - } + bool operator!=(ConstRowIterator const& other) const; /*! * Retrieves an iterator that points to the beginning of the current row. * * @return An iterator that points to the beginning of the current row. */ - ConstIterator begin() const { - return ConstIterator(startValuePtr + *rowPtr, startColumnPtr + *rowPtr); - } + ConstIterator begin() const; /*! * Retrieves an iterator that points past the end of the current row. * * @return An iterator that points past the end of the current row. */ - ConstIterator end() const { - return ConstIterator(startValuePtr + *(rowPtr + 1), startColumnPtr + *(rowPtr + 1)); - } + ConstIterator end() const; private: // The pointer to the value of the first non-zero element of the matrix. @@ -299,74 +264,29 @@ public: * * @param rows The number of rows of the matrix */ - SparseMatrix(uint_fast64_t rows, uint_fast64_t cols) : rowCount(rows), colCount(cols), - nonZeroEntryCount(0), internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { - // Intentionally left empty. - } + SparseMatrix(uint_fast64_t rows, uint_fast64_t cols); /*! * Constructs a square sparse matrix object with the given number rows. * * @param size The number of rows and columns of the matrix. */ - SparseMatrix(uint_fast64_t size = 0) - : rowCount(size), colCount(size), nonZeroEntryCount(0), - internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { - // Intentionally left empty. - } + SparseMatrix(uint_fast64_t size = 0); /*! * Move Constructor. * * @param other The Matrix from which to move the content */ - SparseMatrix(SparseMatrix&& other) - : rowCount(other.rowCount), colCount(other.colCount), nonZeroEntryCount(other.nonZeroEntryCount), - valueStorage(std::move(other.valueStorage)), columnIndications(std::move(other.columnIndications)), - rowIndications(std::move(other.rowIndications)), internalStatus(other.internalStatus), - currentSize(other.currentSize), lastRow(other.lastRow) { - // Now update the source matrix - other.rowCount = 0; - other.colCount = 0; - other.nonZeroEntryCount = 0; - other.internalStatus = MatrixStatus::Error; - other.currentSize = 0; - other.lastRow = 0; - } + SparseMatrix(SparseMatrix&& other); /*! * Copy Constructor. * * @param other The Matrix from which to copy the content */ - SparseMatrix(const SparseMatrix & other) - : rowCount(other.rowCount), colCount(other.colCount), nonZeroEntryCount(other.nonZeroEntryCount), - valueStorage(other.valueStorage), columnIndications(other.columnIndications), - rowIndications(other.rowIndications), internalStatus(other.internalStatus), - currentSize(other.currentSize), lastRow(other.lastRow) { - } - - /*! - * Copy Assignment Operator. - * - * @param other The Matrix from which to copy the content - */ - storm::storage::SparseMatrix& operator=(SparseMatrix const& other) { - this->rowCount = other.rowCount; - this->colCount = other.colCount; - this->nonZeroEntryCount = other.nonZeroEntryCount; + SparseMatrix(const SparseMatrix & other); - this->valueStorage = other.valueStorage; - this->columnIndications = other.columnIndications; - this->rowIndications = other.rowIndications; - - this->internalStatus = other.internalStatus; - this->currentSize = other.currentSize; - this->lastRow = other.lastRow; - - return *this; - } - /*! * Constructs a sparse matrix object with the given (moved) contents. * @@ -379,13 +299,14 @@ public: */ SparseMatrix(uint_fast64_t rowCount, uint_fast64_t colCount, uint_fast64_t nonZeroEntryCount, std::vector&& rowIndications, - std::vector&& columnIndications, std::vector&& values) - : rowCount(rowCount), colCount(colCount), nonZeroEntryCount(nonZeroEntryCount), - valueStorage(values), columnIndications(columnIndications), - rowIndications(rowIndications), internalStatus(MatrixStatus::Initialized), - currentSize(0), lastRow(0) { - // Intentionally left empty. - } + std::vector&& columnIndications, std::vector&& values); + + /*! + * Copy Assignment Operator. + * + * @param other The Matrix from which to copy the content + */ + storm::storage::SparseMatrix& operator=(SparseMatrix const& other); /*! * Initializes the sparse matrix with the given number of non-zero entries @@ -395,30 +316,7 @@ public: * * @param nonZeroEntries The number of non-zero entries this matrix is going to hold. */ - void initialize(uint_fast64_t nonZeroEntries = 0) { - // Check whether initializing the matrix is safe. - if (internalStatus != MatrixStatus::UnInitialized) { - triggerErrorState(); - LOG4CPLUS_ERROR(logger, "Trying to initialize matrix that is not uninitialized."); - throw storm::exceptions::InvalidStateException("Trying to initialize matrix that is not uninitialized."); - } else if ((rowCount * colCount) < nonZeroEntries) { - triggerErrorState(); - LOG4CPLUS_ERROR(logger, "Trying to initialize a matrix with more non-zero entries than there can be."); - throw storm::exceptions::InvalidArgumentException("Trying to initialize a matrix with more non-zero entries than there can be."); - } else { - // If it is safe, initialize necessary members and prepare the - // internal storage. - nonZeroEntryCount = nonZeroEntries; - lastRow = 0; - - if (!prepareInternalStorage()) { - triggerErrorState(); - throw std::bad_alloc(); - } else { - setState(MatrixStatus::Initialized); - } - } - } + void initialize(uint_fast64_t nonZeroEntries = 0); /*! * Sets the matrix element at the given row and column to the given value. After all elements have been added, @@ -433,30 +331,7 @@ public: * @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, T const& value) { - // Check whether the given row and column positions are valid and throw - // error otherwise. - if ((row > rowCount) || (col > colCount)) { - triggerErrorState(); - LOG4CPLUS_ERROR(logger, "Trying to add a value at illegal position (" << row << ", " << col << ") in matrix of size (" << rowCount << ", " << colCount << ")."); - throw storm::exceptions::OutOfRangeException() << "Trying to add a value at illegal position (" << row << ", " << col << ") in matrix of size (" << rowCount << ", " << colCount << ")."; - } - - // 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; - } - lastRow = row; - } - - // Finally, set the element and increase the current size. - valueStorage[currentSize] = value; - columnIndications[currentSize] = col; - - ++currentSize; - } + void addNextValue(const uint_fast64_t row, const uint_fast64_t col, T const& value); /*! * Inserts a value at the given row and column with the given value. After all elements have been inserted, @@ -472,37 +347,7 @@ public: * @param pushRowIndication If set to true, the next row indication value is pushed, otherwise it is added. If the * number of rows was not set in the beginning, then this needs to be true and false otherwise. */ - void insertNextValue(const uint_fast64_t row, const uint_fast64_t col, T const& value, bool pushRowIndication = false) { - // Check whether the given row and column positions are valid and throw - // error otherwise. - if (row < lastRow) { - triggerErrorState(); - LOG4CPLUS_ERROR(logger, "Trying to insert a value at illegal position (" << row << ", " << col << ")."); - throw storm::exceptions::OutOfRangeException() << "Trying to insert a value at illegal position (" << row << ", " << col << ")."; - } - - // If we switched to another row, we have to adjust the missing entries in the rowIndications array. - if (row != lastRow) { - for (uint_fast64_t i = lastRow + 1; i <= row; ++i) { - if (pushRowIndication) { - rowIndications.push_back(currentSize); - } else { - rowIndications[i] = currentSize; - } - } - rowCount = row + 1; - lastRow = row; - } - - // Finally, set the element and increase the current size. - valueStorage.push_back(value); - columnIndications.push_back(col); - ++nonZeroEntryCount; - ++currentSize; - - // Check that we also have the correct number of columns. - colCount = std::max(colCount, col + 1); - } + void insertNextValue(const uint_fast64_t row, const uint_fast64_t col, T const& value, bool pushRowIndication = false); /*! * Inserts an empty row in the matrix. @@ -510,16 +355,7 @@ public: * @param pushRowIndication If set to true, the next row indication value is pushed, otherwise it is added. If the * number of rows was not set in the beginning, then this needs to be true and false otherwise. */ - void insertEmptyRow(bool pushRowIndication = false) { - if (pushRowIndication) { - rowIndications.push_back(currentSize); - } else { - rowIndications[lastRow + 1] = currentSize; - } - - ++rowCount; - ++lastRow; - } + void insertEmptyRow(bool pushRowIndication = false); /* * Finalizes the sparse matrix to indicate that initialization has been completed and the matrix may now be used. @@ -528,35 +364,7 @@ public: * at a fixed location. If the elements have been added to the matrix via insertNextValue, this needs to be true * and false otherwise. */ - void finalize(bool pushSentinelElement = false) { - // Check whether it's safe to finalize the matrix and throw error otherwise. - if (!isInitialized()) { - triggerErrorState(); - LOG4CPLUS_ERROR(logger, "Trying to finalize an uninitialized matrix."); - throw storm::exceptions::InvalidStateException("Trying to finalize an uninitialized matrix."); - } else if (currentSize != nonZeroEntryCount) { - triggerErrorState(); - LOG4CPLUS_ERROR(logger, "Trying to finalize a matrix that was initialized with more non-zero entries than given (expected " << nonZeroEntryCount << " but got " << currentSize << " instead)"); - throw storm::exceptions::InvalidStateException() << "Trying to finalize a matrix that was initialized with more non-zero entries than given (expected " << nonZeroEntryCount << " but got " << currentSize << " instead)."; - } else { - // Fill in the missing entries in the row_indications array. - // (Can happen because of empty rows at the end.) - for (uint_fast64_t i = lastRow + 1; i < rowCount; ++i) { - rowIndications[i] = currentSize; - } - - // 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 rowIndications[i] and rowIndications[i + 1], also for the - // first and last row. - if (pushSentinelElement) { - rowIndications.push_back(nonZeroEntryCount); - } else { - rowIndications[rowCount] = nonZeroEntryCount; - } - - setState(MatrixStatus::ReadReady); - } - } + void finalize(bool pushSentinelElement = false); /*! * Gets the matrix element at the given row and column to the given value. @@ -569,39 +377,7 @@ public: * @returns 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) const { - // Check for illegal access indices. - if ((row > rowCount) || (col > colCount)) { - LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ")."); - throw storm::exceptions::OutOfRangeException("Trying to read a value from illegal position."); - return false; - } - - // 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 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; - } - ++rowStart; - } - - // Set 0 as the content and return false in case the element was not found. - *target = 0; - return false; - } + inline bool getValue(uint_fast64_t row, uint_fast64_t col, T* const target) const; /*! * Gets the matrix element at the given row and column in the form of a reference to it. @@ -615,52 +391,21 @@ public: * * @return A reference to the value at the given position. */ - inline T& getValue(uint_fast64_t row, uint_fast64_t col) { - // Check for illegal access indices. - if ((row > rowCount) || (col > colCount)) { - LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ")."); - throw storm::exceptions::OutOfRangeException("Trying to read a value from illegal position."); - } - - // 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, return it. - if (columnIndications[rowStart] == col) { - return valueStorage[rowStart]; - } - - // 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; - } - ++rowStart; - } - - throw storm::exceptions::InvalidArgumentException("Trying to get a reference to a non-existant value."); - } + inline T& getValue(uint_fast64_t row, uint_fast64_t col); /*! * Returns the number of rows of the matrix. * * @returns The number of rows of the matrix. */ - uint_fast64_t getRowCount() const { - return rowCount; - } + uint_fast64_t getRowCount() const; /*! * Returns the number of columns of the matrix. * * @returns The number of columns of the matrix. */ - uint_fast64_t getColumnCount() const { - return colCount; - } + uint_fast64_t getColumnCount() const; /*! * Checks whether the internal status of the matrix makes it ready for @@ -669,9 +414,7 @@ public: * @returns True iff the internal status of the matrix makes it ready for * reading access. */ - bool isReadReady() { - return (internalStatus == MatrixStatus::ReadReady); - } + bool isReadReady(); /*! * Checks whether the matrix was initialized previously. The matrix may @@ -679,36 +422,28 @@ public: * * @returns True iff the matrix was initialized previously. */ - bool isInitialized() { - return (internalStatus == MatrixStatus::Initialized || internalStatus == MatrixStatus::ReadReady); - } + bool isInitialized(); /*! * Returns the internal state of the matrix. * * @returns The internal state of the matrix. */ - MatrixStatus getState() { - return internalStatus; - } + MatrixStatus getState(); /*! * Checks whether the internal state of the matrix signals an error. * * @returns True iff the internal state of the matrix signals an error. */ - bool hasError() const { - return (internalStatus == MatrixStatus::Error); - } + bool hasError() const; /*! * Returns the number of non-zero entries in the matrix. * * @returns The number of non-zero entries in the matrix. */ - uint_fast64_t getNonZeroEntryCount() const { - return nonZeroEntryCount; - } + uint_fast64_t getNonZeroEntryCount() const; /*! * This function makes the rows given by the bit vector absorbing. @@ -716,13 +451,7 @@ public: * @param rows A bit vector indicating which rows to make absorbing. * @returns True iff the operation was successful. */ - bool makeRowsAbsorbing(storm::storage::BitVector const& rows) { - bool result = true; - for (auto row : rows) { - result &= makeRowAbsorbing(row, row); - } - return result; - } + bool makeRowsAbsorbing(storm::storage::BitVector const& rows); /*! * This function makes the groups of rows given by the bit vector absorbing. @@ -731,16 +460,7 @@ public: * @param rowGroupIndices A vector indicating which rows belong to a given row group. * @return True iff the operation was successful. */ - bool makeRowsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices) { - bool result = true; - for (auto rowGroup : rowGroupConstraint) { - for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) { - result &= makeRowAbsorbing(row, rowGroup); - } - } - return result; - } - + bool makeRowsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices); /*! * This function makes the given row absorbing. This means that all entries will be set to 0 @@ -750,38 +470,7 @@ public: * @param column The index of the column whose value is to be set to 1. * @returns True iff the operation was successful. */ - bool makeRowAbsorbing(const uint_fast64_t row, const uint_fast64_t column) { - // Check whether the accessed state exists. - if (row > rowCount) { - LOG4CPLUS_ERROR(logger, "Trying to make an illegal row " << row << " absorbing."); - throw storm::exceptions::OutOfRangeException() << "Trying to make an illegal row " << row << " absorbing."; - return false; - } - - // Iterate over the elements in the row that are not on the diagonal - // and set them to zero. - uint_fast64_t rowStart = rowIndications[row]; - uint_fast64_t rowEnd = rowIndications[row + 1]; - - // If the row has no elements in it, we cannot make it absorbing, because we would need to - // move all elements in the vector of nonzeros otherwise. - if (rowStart >= rowEnd) { - LOG4CPLUS_ERROR(logger, "Cannot make row " << row << " absorbing, because there is no entry in this row."); - throw storm::exceptions::InvalidStateException() << "Cannot make row " << row << " absorbing, because there is no entry in this row."; - } - - // If there is at least one nonzero entry in this row, we can just set it to one, modify its - // column indication to the one given by the parameter and set all subsequent elements of this - // row to zero. - valueStorage[rowStart] = storm::utility::constGetOne(); - columnIndications[rowStart] = column; - for (uint_fast64_t index = rowStart + 1; index < rowEnd; ++index) { - valueStorage[index] = storm::utility::constGetZero(); - columnIndications[index] = 0; - } - - return true; - } + bool makeRowAbsorbing(const uint_fast64_t row, const uint_fast64_t column); /* * Computes the sum of the elements in the given row whose column bits are set to one on the @@ -792,15 +481,7 @@ public: * @returns The sum of the elements in the given row whose column bits * are set to one on the given constraint. */ - T getConstrainedRowSum(uint_fast64_t row, storm::storage::BitVector const& constraint) const { - T result(0); - for (uint_fast64_t i = rowIndications[row]; i < rowIndications[row + 1]; ++i) { - if (constraint.get(columnIndications[i])) { - result += valueStorage[i]; - } - } - return result; - } + T getConstrainedRowSum(uint_fast64_t row, storm::storage::BitVector const& constraint) const; /*! * Computes a vector whose ith element is the sum of the elements in the ith row where only @@ -814,14 +495,7 @@ public: * those elements are added whose bits are set to true in the given column constraint and only * those rows are treated whose bits are set to true in the given row constraint. */ - std::vector getConstrainedRowSumVector(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint) const { - std::vector result(rowConstraint.getNumberOfSetBits()); - uint_fast64_t currentRowCount = 0; - for (auto row : rowConstraint) { - result[currentRowCount++] = getConstrainedRowSum(row, columnConstraint); - } - return result; - } + std::vector getConstrainedRowSumVector(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint) const; /*! * Computes a vector whose elements represent the sums of selected (given by the column @@ -832,16 +506,7 @@ public: * @param columnConstraint A bit vector that indicates which columns to add. * @returns */ - std::vector getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices, storm::storage::BitVector const& columnConstraint, uint_fast64_t numberOfRows) const { - std::vector result(numberOfRows); - uint_fast64_t currentRowCount = 0; - for (auto rowGroup : rowGroupConstraint) { - for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) { - result[currentRowCount++] = getConstrainedRowSum(row, columnConstraint); - } - } - return result; - } + std::vector getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices, storm::storage::BitVector const& columnConstraint, uint_fast64_t numberOfRows) const; /*! * Creates a submatrix of the current matrix by dropping all rows and columns whose bits are not @@ -851,63 +516,7 @@ public: * @returns A matrix corresponding to a submatrix of the current matrix in which only rows and * columns given by the constraint are kept and all others are dropped. */ - SparseMatrix getSubmatrix(storm::storage::BitVector const& constraint) const { - LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix with " << constraint.getNumberOfSetBits() << " rows."); - - // Check for valid constraint. - if (constraint.getNumberOfSetBits() == 0) { - LOG4CPLUS_ERROR(logger, "Trying to create a sub-matrix of size 0."); - throw storm::exceptions::InvalidArgumentException("Trying to create a sub-matrix of size 0."); - } - - // First, we need to determine the number of non-zero entries of the - // sub-matrix. - uint_fast64_t subNonZeroEntries = 0; - for (auto rowIndex : constraint) { - for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) { - if (constraint.get(columnIndications[i])) { - ++subNonZeroEntries; - } - } - } - - // Create and initialize resulting matrix. - SparseMatrix result(constraint.getNumberOfSetBits()); - result.initialize(subNonZeroEntries); - - // Create a temporary vecotr that stores for each index whose bit is set - // to true the number of bits that were set before that particular index. - std::vector bitsSetBeforeIndex; - bitsSetBeforeIndex.reserve(colCount); - - // Compute the information to fill this vector. - uint_fast64_t lastIndex = 0; - uint_fast64_t currentNumberOfSetBits = 0; - for (auto index : constraint) { - while (lastIndex <= index) { - bitsSetBeforeIndex.push_back(currentNumberOfSetBits); - ++lastIndex; - } - ++currentNumberOfSetBits; - } - - // Copy over selected entries. - uint_fast64_t rowCount = 0; - for (auto rowIndex : constraint) { - for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) { - if (constraint.get(columnIndications[i])) { - result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[i]], valueStorage[i]); - } - } - - ++rowCount; - } - - // Finalize sub-matrix and return result. - result.finalize(); - LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); - return result; - } + SparseMatrix getSubmatrix(storm::storage::BitVector const& constraint) const; /*! * Creates a submatrix of the current matrix by keeping only row groups and columns in the given @@ -918,217 +527,34 @@ public: * @returns A matrix corresponding to a submatrix of the current matrix in which only row groups * and columns given by the row group constraint are kept and all others are dropped. */ - SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices) const { - LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size)."); - - // First, we need to determine the number of non-zero entries and the number of rows of the sub-matrix. - uint_fast64_t subNonZeroEntries = 0; - uint_fast64_t subRowCount = 0; - for (auto index : rowGroupConstraint) { - subRowCount += rowGroupIndices[index + 1] - rowGroupIndices[index]; - for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { - for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { - if (rowGroupConstraint.get(columnIndications[j])) { - ++subNonZeroEntries; - } - } - } - } - - LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << rowGroupConstraint.getNumberOfSetBits() << "."); - - // Create and initialize resulting matrix. - SparseMatrix result(subRowCount, rowGroupConstraint.getNumberOfSetBits()); - result.initialize(subNonZeroEntries); - - // Create a temporary vector that stores for each index whose bit is set - // to true the number of bits that were set before that particular index. - std::vector bitsSetBeforeIndex; - bitsSetBeforeIndex.reserve(colCount); - - // Compute the information to fill this vector. - uint_fast64_t lastIndex = 0; - uint_fast64_t currentNumberOfSetBits = 0; - for (auto index : rowGroupConstraint) { - while (lastIndex <= index) { - bitsSetBeforeIndex.push_back(currentNumberOfSetBits); - ++lastIndex; - } - ++currentNumberOfSetBits; - } - - // Copy over selected entries. - uint_fast64_t rowCount = 0; - for (auto index : rowGroupConstraint) { - for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { - for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { - if (rowGroupConstraint.get(columnIndications[j])) { - result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]); - } - } - ++rowCount; - } - } - - // Finalize sub-matrix and return result. - result.finalize(); - LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); - return result; - } + SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices) const; - SparseMatrix getSubmatrix(std::vector const& rowGroupToRowIndexMapping, std::vector const& rowGroupIndices, bool insertDiagonalEntries = true) const { - LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size)."); - - // First, we need to count how many non-zero entries the resulting matrix will have and reserve space for diagonal - // entries. - uint_fast64_t subNonZeroEntries = 0; - for (uint_fast64_t rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) { - // Determine which row we need to select from the current row group. - uint_fast64_t rowToCopy = rowGroupIndices[rowGroupIndex] + rowGroupToRowIndexMapping[rowGroupIndex]; - - // Iterate through that row and count the number of slots we have to reserve for copying. - bool foundDiagonalElement = false; - for (uint_fast64_t i = rowIndications[rowToCopy], rowEnd = rowIndications[rowToCopy + 1]; i < rowEnd; ++i) { - if (columnIndications[i] == rowGroupIndex) { - foundDiagonalElement = true; - } - ++subNonZeroEntries; - } - if (insertDiagonalEntries && !foundDiagonalElement) { - ++subNonZeroEntries; - } - } - - LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << (rowGroupIndices.size() - 1) << "x" << colCount << " with " << subNonZeroEntries << " non-zero elements."); - - // Now create the matrix to be returned with the appropriate size. - SparseMatrix submatrix(rowGroupIndices.size() - 1, colCount); - submatrix.initialize(subNonZeroEntries); - - // Copy over the selected lines from the source matrix. - for (uint_fast64_t rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) { - // Determine which row we need to select from the current row group. - uint_fast64_t rowToCopy = rowGroupIndices[rowGroupIndex] + rowGroupToRowIndexMapping[rowGroupIndex]; - - // Iterate through that row and copy the entries. This also inserts a zero element on the diagonal if there - // is no entry yet. - bool insertedDiagonalElement = false; - for (uint_fast64_t i = rowIndications[rowToCopy], rowEnd = rowIndications[rowToCopy + 1]; i < rowEnd; ++i) { - if (columnIndications[i] == rowGroupIndex) { - insertedDiagonalElement = true; - } else if (insertDiagonalEntries && !insertedDiagonalElement && columnIndications[i] > rowGroupIndex) { - submatrix.addNextValue(rowGroupIndex, rowGroupIndex, storm::utility::constGetZero()); - insertedDiagonalElement = true; - } - submatrix.addNextValue(rowGroupIndex, columnIndications[i], valueStorage[i]); - } - if (insertDiagonalEntries && !insertedDiagonalElement) { - submatrix.addNextValue(rowGroupIndex, rowGroupIndex, storm::utility::constGetZero()); - } - } - - // Finalize created matrix and return result. - submatrix.finalize(); - LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); - return submatrix; - } + SparseMatrix getSubmatrix(std::vector const& rowGroupToRowIndexMapping, std::vector const& rowGroupIndices, bool insertDiagonalEntries = true) const; /*! * Performs a change to the matrix that is needed if this matrix is to be interpreted as a * set of linear equations. In particular, it transforms A to (1-A), meaning that the elements * on the diagonal are inverted with respect to addition and the other elements are negated. */ - void convertToEquationSystem() { - invertDiagonal(); - negateAllNonDiagonalElements(); - } + void convertToEquationSystem(); /*! * Inverts all elements on the diagonal, i.e. sets the diagonal values to 1 minus their previous * value. Requires the matrix to contain each diagonal element and to be square. */ - void invertDiagonal() { - // Check if the matrix is square, because only then it makes sense to perform this - // transformation. - if (this->getRowCount() != this->getColumnCount()) { - throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to be square!"; - } - - // Now iterate over all rows and set the diagonal elements to the inverted value. - // If there is a row without the diagonal element, an exception is thrown. - T one = storm::utility::constGetOne(); - bool foundDiagonalElement = false; - for (uint_fast64_t row = 0; row < rowCount; ++row) { - uint_fast64_t rowStart = rowIndications[row]; - uint_fast64_t rowEnd = rowIndications[row + 1]; - foundDiagonalElement = false; - while (rowStart < rowEnd) { - if (columnIndications[rowStart] == row) { - valueStorage[rowStart] = one - valueStorage[rowStart]; - foundDiagonalElement = true; - break; - } - ++rowStart; - } - - // Throw an exception if a row did not have an element on the diagonal. - if (!foundDiagonalElement) { - throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to contain all diagonal entries!"; - } - } - } + void invertDiagonal(); /*! * Negates all non-zero elements that are not on the diagonal. */ - void negateAllNonDiagonalElements() { - // Check if the matrix is square, because only then it makes sense to perform this - // transformation. - if (this->getRowCount() != this->getColumnCount()) { - throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::invertDiagonal requires the Matrix to be square!"; - } - - // Iterate over all rows and negate all the elements that are not on the diagonal. - for (uint_fast64_t row = 0; row < rowCount; ++row) { - uint_fast64_t rowStart = rowIndications[row]; - uint_fast64_t rowEnd = rowIndications[row + 1]; - while (rowStart < rowEnd) { - if (columnIndications[rowStart] != row) { - valueStorage[rowStart] = - valueStorage[rowStart]; - } - ++rowStart; - } - } - } + void negateAllNonDiagonalElements(); /*! * Calculates the Jacobi-Decomposition of this sparse matrix. * The source Sparse Matrix must be square. * @return A std::pair containing the matrix L+U and the inverted diagonal matrix D^-1 */ - SparseJacobiDecomposition_t getJacobiDecomposition() const { - uint_fast64_t rowCount = this->getRowCount(); - uint_fast64_t colCount = this->getColumnCount(); - if (rowCount != colCount) { - throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::getJacobiDecomposition requires the Matrix to be square."; - } - storm::storage::SparseMatrix resultLU(*this); - storm::storage::SparseMatrix resultDinv(rowCount, colCount); - // no entries apart from those on the diagonal (rowCount) - resultDinv.initialize(rowCount); - - // constant 1 for diagonal inversion - T constOne = storm::utility::constGetOne(); - - // copy diagonal entries to other matrix - for (uint_fast64_t i = 0; i < rowCount; ++i) { - resultDinv.addNextValue(i, i, constOne / resultLU.getValue(i, i)); - resultLU.getValue(i, i) = storm::utility::constGetZero(); - } - resultDinv.finalize(); - - return std::make_pair(std::move(resultLU), std::move(resultDinv)); - } + SparseJacobiDecomposition_t getJacobiDecomposition() const; /*! * Performs a pointwise matrix multiplication of the matrix with the given matrix and returns a @@ -1140,29 +566,7 @@ public: * @returns A vector containing the sum of the elements in each row of the matrix resulting from * pointwise multiplication of the current matrix with the given matrix. */ - std::vector getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const { - // Prepare result. - std::vector result(rowCount, storm::utility::constGetZero()); - - // Iterate over all elements of the current matrix and either continue with the next element - // in case the given matrix does not have a non-zero element at this column position, or - // multiply the two entries and add the result to the corresponding position in the vector. - for (uint_fast64_t row = 0; row < rowCount && row < otherMatrix.rowCount; ++row) { - for (uint_fast64_t element = rowIndications[row], nextOtherElement = otherMatrix.rowIndications[row]; element < rowIndications[row + 1] && nextOtherElement < otherMatrix.rowIndications[row + 1]; ++element) { - if (columnIndications[element] < otherMatrix.columnIndications[nextOtherElement]) { - continue; - } else { - // If the precondition of this method (i.e. that the given matrix is a submatrix - // of the current one) was fulfilled, we know now that the two elements are in - // the same column, so we can multiply and add them to the row sum vector. - result[row] += otherMatrix.valueStorage[nextOtherElement] * valueStorage[element]; - ++nextOtherElement; - } - } - } - - return result; - } + std::vector getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; /*! * Multiplies the matrix with the given vector and writes the result to given result vector. @@ -1173,37 +577,14 @@ public: * @returns The product of the matrix and the given vector as the content of the given result * vector. */ - void multiplyWithVector(std::vector const& vector, std::vector& result) const { -#ifdef STORM_HAVE_INTELTBB - tbb::parallel_for(tbb::blocked_range(0, result.size()), tbbHelper_MatrixRowVectorScalarProduct, std::vector, T>(this, &vector, &result)); -#else - ConstRowIterator rowIt = this->begin(); - - for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++rowIt) { - *resultIt = storm::utility::constGetZero(); - - for (auto elementIt = rowIt.begin(), elementIte = rowIt.end(); elementIt != elementIte; ++elementIt) { - *resultIt += elementIt.value() * vector[elementIt.column()]; - } - } -#endif - } + void multiplyWithVector(std::vector const& vector, std::vector& result) const; /*! * Returns the size of the matrix in memory measured in bytes. * * @returns The size of the matrix in memory measured in bytes. */ - uint_fast64_t getSizeInMemory() const { - uint_fast64_t size = sizeof(*this); - // Add value_storage size. - size += sizeof(T) * valueStorage.capacity(); - // Add column_indications size. - size += sizeof(uint_fast64_t) * columnIndications.capacity(); - // Add row_indications size. - size += sizeof(uint_fast64_t) * rowIndications.capacity(); - return size; - } + uint_fast64_t getSizeInMemory() const; /*! * Returns an object representing the consecutive rows given by the parameters. @@ -1212,9 +593,7 @@ public: * @param endRow The ending row (which is included in the result). * @return An object representing the consecutive rows given by the parameters. */ - Rows getRows(uint_fast64_t startRow, uint_fast64_t endRow) const { - return Rows(this->valueStorage.data() + this->rowIndications[startRow], this->columnIndications.data() + this->rowIndications[startRow], this->rowIndications[endRow + 1] - this->rowIndications[startRow]); - } + Rows getRows(uint_fast64_t startRow, uint_fast64_t endRow) const; /*! * Returns an object representing the given row. @@ -1222,9 +601,7 @@ public: * @param row The chosen row. * @return An object representing the given row. */ - Rows getRow(uint_fast64_t row) const { - return getRows(row, row); - } + Rows getRow(uint_fast64_t row) const; /*! * Returns a const iterator to the rows of the matrix. @@ -1232,18 +609,14 @@ public: * @param initialRow The initial row to which this iterator points. * @return A const iterator to the rows of the matrix. */ - ConstRowIterator begin(uint_fast64_t initialRow = 0) const { - return ConstRowIterator(this->valueStorage.data(), this->columnIndications.data(), this->rowIndications.data() + initialRow); - } + ConstRowIterator begin(uint_fast64_t initialRow = 0) const; /*! * Returns a const iterator that points to past the last row of the matrix. * * @return A const iterator that points to past the last row of the matrix */ - ConstRowIterator end() const { - return ConstRowIterator(this->valueStorage.data(), this->columnIndications.data(), this->rowIndications.data() + rowCount); - } + ConstRowIterator end() const; /*! * Returns a const iterator that points past the given row. @@ -1251,9 +624,7 @@ public: * @param row The row past which this iterator points. * @return A const iterator that points past the given row. */ - ConstRowIterator end(uint_fast64_t row) const { - return ConstRowIterator(this->valueStorage.data(), this->columnIndications.data(), this->rowIndications.data() + row + 1); - } + ConstRowIterator end(uint_fast64_t row) const; /*! * Returns an iterator to the columns of the non-zero entries of the matrix. @@ -1261,18 +632,14 @@ public: * @param row If given, the iterator will start at the specified row. * @returns An iterator to the columns of the non-zero entries of the matrix. */ - ConstIndexIterator constColumnIteratorBegin(uint_fast64_t row = 0) const { - return &(this->columnIndications[0]) + this->rowIndications[row]; - } + ConstIndexIterator constColumnIteratorBegin(uint_fast64_t row = 0) const; /*! * Returns an iterator that points to the first element after the matrix. * * @returns An iterator that points to the first element after the matrix. */ - ConstIndexIterator constColumnIteratorEnd() const { - return &(this->columnIndications[0]) + this->rowIndications[rowCount]; - } + ConstIndexIterator constColumnIteratorEnd() const; /*! * Returns an iterator that points to the first element after the given row. @@ -1280,9 +647,7 @@ public: * @param row The row past which this iterator has to point. * @returns An iterator that points to the first element after the matrix. */ - ConstIndexIterator constColumnIteratorEnd(uint_fast64_t row) const { - return &(this->columnIndications[0]) + this->rowIndications[row + 1]; - } + ConstIndexIterator constColumnIteratorEnd(uint_fast64_t row) const; /*! * Returns an iterator to the values of the non-zero entries of the matrix. @@ -1290,18 +655,14 @@ public: * @param row If given, the iterator will start at the specified row. * @returns An iterator to the values of the non-zero entries of the matrix. */ - ConstValueIterator constValueIteratorBegin(uint_fast64_t row = 0) const { - return &(this->valueStorage[0]) + this->rowIndications[row]; - } + ConstValueIterator constValueIteratorBegin(uint_fast64_t row = 0) const; /*! * Returns an iterator that points to the first element after the matrix. * * @returns An iterator that points to the first element after the matrix. */ - ConstValueIterator constValueIteratorEnd() const { - return &(this->valueStorage[0]) + this->rowIndications[rowCount]; - } + ConstValueIterator constValueIteratorEnd() const; /*! * Returns an iterator that points to the first element after the given row. @@ -1309,9 +670,7 @@ public: * @param row The row past which this iterator has to point. * @returns An iterator that points to the first element after the matrix. */ - ConstValueIterator constValueIteratorEnd(uint_fast64_t row) const { - return &(this->valueStorage[0]) + this->rowIndications[row + 1]; - } + ConstValueIterator constValueIteratorEnd(uint_fast64_t row) const; /*! * Computes the sum of the elements in a given row. @@ -1319,13 +678,7 @@ public: * @param row The row that should be summed. * @return Sum of the row. */ - T getRowSum(uint_fast64_t row) const { - T sum = storm::utility::constGetZero(); - for (auto it = this->constValueIteratorBegin(row), ite = this->constValueIteratorEnd(row); it != ite; ++it) { - sum += *it; - } - return sum; - } + T getRowSum(uint_fast64_t row) const; /*! * Checks if the current matrix is a submatrix of the given matrix, where a matrix A is called a @@ -1335,39 +688,14 @@ public: * @param matrix The matrix that possibly is a "supermatrix" of the current matrix. * @returns True iff the current matrix is a submatrix of the given matrix. */ - bool isSubmatrixOf(SparseMatrix const& matrix) const { - // Check for matching sizes. - if (this->getRowCount() != matrix.getRowCount()) return false; - if (this->getColumnCount() != matrix.getColumnCount()) return false; - - // Check the subset property for all rows individually. - for (uint_fast64_t row = 0; row < this->getRowCount(); ++row) { - for (uint_fast64_t elem = rowIndications[row], elem2 = matrix.rowIndications[row]; elem < rowIndications[row + 1] && elem < matrix.rowIndications[row + 1]; ++elem) { - // Skip over all entries of the other matrix that are before the current entry in - // the current matrix. - while (elem2 < matrix.rowIndications[row + 1] && matrix.columnIndications[elem2] < columnIndications[elem]) { - ++elem2; - } - if (!(elem2 < matrix.rowIndications[row + 1]) || columnIndications[elem] != matrix.columnIndications[elem2]) { - return false; - } - } - } - return true; - } + bool isSubmatrixOf(SparseMatrix const& matrix) const; /*! * Retrieves a compressed string representation of the matrix. * * @returns a compressed string representation of the matrix. */ - std::string toStringCompressed() const { - std::stringstream result; - result << rowIndications << std::endl; - result << columnIndications << std::endl; - result << valueStorage << std::endl; - return result.str(); - } + std::string toStringCompressed() const; /*! * Retrieves a (non-compressed) string representation of the matrix. @@ -1378,82 +706,13 @@ public: * given, rows of different groups will be separated by a dashed line. * @returns A (non-compressed) string representation of the matrix. */ - std::string toString(std::vector const* rowGroupIndices = nullptr) const { - std::stringstream result; - uint_fast64_t currentNondeterministicChoiceIndex = 0; - - // Print column numbers in header. - result << "\t\t"; - for (uint_fast64_t i = 0; i < colCount; ++i) { - result << i << "\t"; - } - result << std::endl; - - // Iterate over all rows. - for (uint_fast64_t i = 0; i < rowCount; ++i) { - uint_fast64_t nextIndex = rowIndications[i]; - - // If we need to group rows, print a dashed line in case we have moved to the next group of rows. - if (rowGroupIndices != nullptr) { - if (i == (*rowGroupIndices)[currentNondeterministicChoiceIndex]) { - if (i != 0) { - result << "\t(\t"; - for (uint_fast64_t j = 0; j < colCount - 2; ++j) { - result << "----"; - if (j == 1) { - result << "\t" << currentNondeterministicChoiceIndex << "\t"; - } - } - result << "\t)" << std::endl; - } - ++currentNondeterministicChoiceIndex; - } - } - - // Print the actual row. - result << i << "\t(\t"; - uint_fast64_t currentRealIndex = 0; - while (currentRealIndex < colCount) { - if (nextIndex < rowIndications[i + 1] && currentRealIndex == columnIndications[nextIndex]) { - result << valueStorage[nextIndex] << "\t"; - ++nextIndex; - } else { - result << "0\t"; - } - ++currentRealIndex; - } - result << "\t)\t" << i << std::endl; - } - - // Print column numbers in footer. - result << "\t\t"; - for (uint_fast64_t i = 0; i < colCount; ++i) { - result << i << "\t"; - } - result << std::endl; - - // Return final result. - return result.str(); - } + std::string toString(std::vector const* rowGroupIndices = nullptr) const; /*! * Calculates a hash over all values contained in this Sparse Matrix. * @return size_t A Hash Value */ - std::size_t getHash() const { - std::size_t result = 0; - - boost::hash_combine(result, rowCount); - boost::hash_combine(result, colCount); - boost::hash_combine(result, nonZeroEntryCount); - boost::hash_combine(result, currentSize); - boost::hash_combine(result, lastRow); - boost::hash_combine(result, storm::utility::Hash::getHash(valueStorage)); - boost::hash_combine(result, storm::utility::Hash::getHash(columnIndications)); - boost::hash_combine(result, storm::utility::Hash::getHash(rowIndications)); - - return result; - } + std::size_t getHash() const; private: @@ -1509,18 +768,14 @@ private: /*! * Sets the internal status to signal an error. */ - void triggerErrorState() { - setState(MatrixStatus::Error); - } + void triggerErrorState(); /*! * 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) { - internalStatus = (internalStatus == MatrixStatus::Error) ? internalStatus : new_state; - } + void setState(const MatrixStatus new_state); /*! * Prepares the internal storage. It relies on the number of non-zero entries and the @@ -1530,29 +785,14 @@ private: * @param initializeElements If set to true, all entries are initialized. * @return True on success, false otherwise (allocation failed). */ - bool prepareInternalStorage(bool initializeElements = true) { - if (initializeElements) { - // Set up the arrays for the elements that are not on the diagonal. - valueStorage.resize(nonZeroEntryCount, storm::utility::constGetZero()); - columnIndications.resize(nonZeroEntryCount, 0); - - // Set up the rowIndications vector and reserve one element more than there are rows in - // order to put a sentinel element at the end, which eases iteration process. - rowIndications.resize(rowCount + 1, 0); - - // Return whether all the allocations could be made without error. - return ((valueStorage.capacity() >= nonZeroEntryCount) && (columnIndications.capacity() >= nonZeroEntryCount) - && (rowIndications.capacity() >= (rowCount + 1))); - } else { - // If it was not requested to initialize the elements, we simply reserve the space. - valueStorage.reserve(nonZeroEntryCount); - columnIndications.reserve(nonZeroEntryCount); - rowIndications.reserve(rowCount + 1); - return true; - } - } + bool prepareInternalStorage(bool initializeElements = true); }; +// Extern template declaration to tell the compiler that there already is an instanciation of this template somewhere. +// The extern instance will be found by the linker. Prevents multiple instantiations of the same template specialisation. +extern template class SparseMatrix; +extern template class SparseMatrix; + #ifdef STORM_HAVE_INTELTBB /*! * This function is a helper for Parallel Execution of the multipliyWithVector functionality. @@ -1560,32 +800,20 @@ private: */ template class tbbHelper_MatrixRowVectorScalarProduct { + public: + tbbHelper_MatrixRowVectorScalarProduct(M const* matrixA, V const* vectorX, V * resultVector); + + void operator() (const tbb::blocked_range& r) const; + private: V * resultVector; V const* vectorX; M const* matrixA; - - public: - tbbHelper_MatrixRowVectorScalarProduct(M const* matrixA, V const* vectorX, V * resultVector) : matrixA(matrixA), vectorX(vectorX), resultVector(resultVector) {} - - void operator() (const tbb::blocked_range& r) const { - for (uint_fast64_t row = r.begin(); row < r.end(); ++row) { - uint_fast64_t index = matrixA->rowIndications.at(row); - uint_fast64_t indexEnd = matrixA->rowIndications.at(row + 1); - - // Initialize the result to be 0. - T element = storm::utility::constGetZero(); - - for (; index != indexEnd; ++index) { - element += matrixA->valueStorage.at(index) * vectorX->at(matrixA->columnIndications.at(index)); - } - - // Write back to the result Vector - resultVector->at(row) = element; - } - } }; + + // Extern template declaration to tell the compiler that there already is an instanciation of this template somewhere. + extern template class tbbHelper_MatrixRowVectorScalarProduct, std::vector, double>; #endif diff --git a/src/storage/StronglyConnectedComponentDecomposition.h b/src/storage/StronglyConnectedComponentDecomposition.h index a2f3e365b..e81221398 100644 --- a/src/storage/StronglyConnectedComponentDecomposition.h +++ b/src/storage/StronglyConnectedComponentDecomposition.h @@ -1,10 +1,11 @@ #ifndef STORM_STORAGE_STRONGLYCONNECTEDCOMPONENTDECOMPOSITION_H_ #define STORM_STORAGE_STRONGLYCONNECTEDCOMPONENTDECOMPOSITION_H_ +#include "src/storage/SparseMatrix.h" #include "src/storage/Decomposition.h" #include "src/storage/VectorSet.h" #include "src/storage/BitVector.h" -#include "src/storage/SparseMatrix.h" + namespace storm { namespace models { diff --git a/src/storm.cpp b/src/storm.cpp index a0d8006c9..6e655c0fb 100644 --- a/src/storm.cpp +++ b/src/storm.cpp @@ -21,9 +21,9 @@ #include #include "storm-config.h" +#include "src/storage/SparseMatrix.h" #include "src/models/Dtmc.h" #include "src/models/MarkovAutomaton.h" -#include "src/storage/SparseMatrix.h" #include "src/models/AtomicPropositionsLabeling.h" #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"