diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index aa8e9f3ea..f2ffb95b3 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -1,11 +1,3 @@ -/* - * SparseMatrix.cpp - * - * Created on: Nov 18, 2013 - * Author: Manuel Sascha Weiand - */ - -#include #include #include "src/storage/SparseMatrix.h" @@ -16,1127 +8,981 @@ extern log4cplus::Logger logger; 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); - } - } - } + namespace storage { + + template + template + SparseMatrix::BaseIterator::BaseIterator(ValueType* valuePtr, uint_fast64_t const* columnPtr) : valuePtr(valuePtr), columnPtr(columnPtr) { + // Intentionally left empty. + } - 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 << ")."; - } + template + template + SparseMatrix::BaseIterator::BaseIterator(SparseMatrix::BaseIterator const& other) : valuePtr(other.valuePtr), columnPtr(other.columnPtr) { + // Intentionally left empty. + } - // 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; - } + template + template + typename SparseMatrix::template BaseIterator& SparseMatrix::BaseIterator::operator=(BaseIterator const& other) { + if (this != &other) { + valuePtr = other.valuePtr, + columnPtr = other.columnPtr; + } + return *this; + } - // Finally, set the element and increase the current size. - valueStorage[currentSize] = value; - columnIndications[currentSize] = col; + template + template + SparseMatrix::BaseIterator& SparseMatrix::BaseIterator::operator++() { + this->valuePtr++; + this->columnPtr++; + return *this; + } + + template + template + SparseMatrix::BaseIterator& SparseMatrix::BaseIterator::operator*() { + return *this; + } + + template + template + bool SparseMatrix::BaseIterator::operator!=(BaseIterator const& other) const { + return this->valuePtr != other.valuePtr; + } + + template + template + bool SparseMatrix::BaseIterator::operator==(BaseIterator const& other) const { + return this->valuePtr == other.valuePtr; + } - ++currentSize; - } + template + template + uint_fast64_t SparseMatrix::BaseIterator::column() const { + return *columnPtr; + } + + template + template + ValueType& SparseMatrix::BaseIterator::value() const { + return *valuePtr; + } + + template + SparseMatrix::rows::rows(T* valuePtr, uint_fast64_t const* columnPtr, uint_fast64_t entryCount) : valuePtr(valuePtr), columnPtr(columnPtr), entryCount(entryCount) { + // Intentionally left empty. + } + + template + typename SparseMatrix::iterator SparseMatrix::rows::begin() { + return iterator(valuePtr, columnPtr); + } + + template + typename SparseMatrix::iterator SparseMatrix::rows::end() { + return iterator(valuePtr + entryCount, columnPtr + entryCount); + } - 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 << ")."; - } + template + SparseMatrix::const_rows::const_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::const_iterator SparseMatrix::const_rows::begin() const { + return const_iterator(valuePtr, columnPtr); + } + + template + typename SparseMatrix::const_iterator SparseMatrix::const_rows::end() const { + return const_iterator(valuePtr + entryCount, columnPtr + entryCount); + } + + template + SparseMatrix::SparseMatrix(uint_fast64_t rows, uint_fast64_t columns, uint_fast64_t entries) : rowCount(rows), columnCount(columns), entryCount(entries), internalStatus(UNINITIALIZED), valueStorage(), columnIndications(), rowIndications(), currentEntryCount(0), lastRow(0), lastColumn(0) { + storagePreallocated = rows != 0 && columns != 0 && entries != 0; + prepareInternalStorage(); + } + + template + SparseMatrix::SparseMatrix(uint_fast64_t size, uint_fast64_t entries) : SparseMatrix(size, size, entries) { + // Intentionally left empty. + } + + template + SparseMatrix::SparseMatrix(SparseMatrix const& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), storagePreallocated(other.storagePreallocated), valueStorage(other.valueStorage), columnIndications(other.columnIndications), rowIndications(other.rowIndications), internalStatus(other.internalStatus), currentEntryCount(other.currentEntryCount), lastRow(other.lastRow), lastColumn(other.lastColumn) { + // Intentionally left empty. + } - // 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); + template + SparseMatrix::SparseMatrix(SparseMatrix&& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), storagePreallocated(other.storagePreallocated), valueStorage(std::move(other.valueStorage)), columnIndications(std::move(other.columnIndications)), rowIndications(std::move(other.rowIndications)), internalStatus(other.internalStatus), currentEntryCount(other.currentEntryCount), lastRow(other.lastRow), lastColumn(other.lastColumn) { + // Now update the source matrix + other.rowCount = 0; + other.columnCount = 0; + other.entryCount = 0; + other.storagePreallocated = false; + other.internalStatus = MatrixStatus::UNINITIALIZED; + other.currentEntryCount = 0; + other.lastRow = 0; + other.lastColumn = 0; + } + + template + SparseMatrix::SparseMatrix(uint_fast64_t columnCount, std::vector&& rowIndications, std::vector&& columnIndications, std::vector&& values) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(values.size()), valueStorage(std::move(values)), columnIndications(std::move(columnIndications)), rowIndications(std::move(rowIndications)), internalStatus(INITIALIZED), currentEntryCount(0), lastRow(0), lastColumn(0) { + // Intentionally left empty. + } + + template + SparseMatrix& SparseMatrix::operator=(SparseMatrix const& other) { + // Only perform assignment if source and target are not the same. + if (this != &other) { + rowCount = other.rowCount; + columnCount = other.columnCount; + entryCount = other.entryCount; + + valueStorage = other.valueStorage; + columnIndications = other.columnIndications; + rowIndications = other.rowIndications; + + internalStatus = other.internalStatus; + currentEntryCount = other.currentEntryCount; + lastRow = other.lastRow; + lastColumn = other.lastColumn; + } + + return *this; + } + + template + void SparseMatrix::addNextValue(uint_fast64_t row, uint_fast64_t column, T const& value) { + // Depending on whether the internal data storage was preallocated or not, adding the value is done somewhat + // differently. + if (storagePreallocated) { + // Check whether the given row and column positions are valid and throw error otherwise. + if (row > rowCount || column > columnCount) { + throw storm::exceptions::OutOfRangeException() << "Illegal call to SparseMatrix::addNextValue: adding entry at out-of-bounds position (" << row << ", " << column << ") in matrix of size (" << rowCount << ", " << columnCount << ")."; + } + } + + // Check that we did not move backwards wrt. the row. + if (row < lastRow) { + throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrix::addNextValue: adding an element in row " << row << ", but an element in row " << lastRow << " has already been added." << std::endl; + } + + // Check that we did not move backwards wrt. to column. + if (row == lastRow && column < lastColumn) { + throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrix::addNextValue: adding an element in column " << column << " in row " << row << ", but an element in column " << lastColumn << " has already been added in that row." << std::endl; + } + + // If we switched to another row, we have to adjust the missing entries in the row indices vector. + if (row != lastRow) { + if (storagePreallocated) { + // If the storage was preallocated, we can access the elements in the vectors with the subscript + // operator. + for (uint_fast64_t i = lastRow + 1; i <= row; ++i) { + rowIndications[i] = currentEntryCount; + } } else { - rowIndications[i] = currentSize; + // Otherwise, we need to push the correct values to the vectors, which might trigger reallocations. + for (uint_fast64_t i = lastRow + 1; i <= row; ++i) { + rowIndications.push_back(currentEntryCount); + } } - } - 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); + lastRow = row; + } + + lastColumn = column; + + // Finally, set the element and increase the current size. + if (storagePreallocated) { + valueStorage[currentEntryCount] = value; + columnIndications[currentEntryCount] = column; + ++currentEntryCount; } else { - rowIndications[rowCount] = nonZeroEntryCount; + valueStorage.push_back(value); + columnIndications.push_back(column); } + } + + template + void SparseMatrix::finalize() { + // Check whether it's safe to finalize the matrix and throw error otherwise. + if (internalStatus == INITIALIZED) { + throw storm::exceptions::InvalidStateException() << "Illegal call to SparseMatrix::finalize: finalizing an initialized matrix is forbidden."; + } else if (storagePreallocated && currentEntryCount != entryCount) { + throw storm::exceptions::InvalidStateException() << "Illegal call to SparseMatrix::finalize: expected " << entryCount << " entries, but got " << currentEntryCount << " instead."; + } else { + // Fill in the missing entries in the row indices array, as there may be empty rows at the end. + if (storagePreallocated) { + for (uint_fast64_t i = lastRow + 1; i < rowCount; ++i) { + rowIndications[i] = currentEntryCount; + } + } else { + for (uint_fast64_t i = lastRow + 1; i < rowCount; ++i) { + rowIndications.push_back(currentEntryCount); + } + } + + // We put a sentinel element at the last position of the row indices 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 (storagePreallocated) { + rowIndications[rowCount] = entryCount; + } else { + rowIndications.push_back(entryCount); + } - 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::constantOne(); - columnIndications[rowStart] = column; - for (uint_fast64_t index = rowStart + 1; index < rowEnd; ++index) { - valueStorage[index] = storm::utility::constantZero(); - 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, bool insertDiagonalEntries) const { - return getSubmatrix(rowGroupConstraint, rowGroupConstraint, rowGroupIndices, insertDiagonalEntries); - } - - template - SparseMatrix SparseMatrix::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector const& rowGroupIndices, bool insertDiagonalEntries) 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) { - bool foundDiagonalElement = false; + internalStatus = INITIALIZED; + } + } + + template + uint_fast64_t SparseMatrix::getRowCount() const { + return rowCount; + } + + template + uint_fast64_t SparseMatrix::getColumnCount() const { + return columnCount; + } + + template + bool SparseMatrix::isInitialized() { + return internalStatus == INITIALIZED; + } + + template + uint_fast64_t SparseMatrix::getEntryCount() const { + return entryCount; + } + + template + void SparseMatrix::makeRowsAbsorbing(storm::storage::BitVector const& rows) { + for (auto row : rows) { + makeRowAbsorbing(row, row); + } + } + + template + void SparseMatrix::makeRowsAbsorbing(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices) { + for (auto rowGroup : rowGroupConstraint) { + for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) { + makeRowAbsorbing(row, rowGroup); + } + } + } + + template + void SparseMatrix::makeRowAbsorbing(const uint_fast64_t row, const uint_fast64_t column) { + if (row > rowCount) { + throw storm::exceptions::OutOfRangeException() << "Illegal call to SparseMatrix::makeRowAbsorbing: access to row " << row << " is out of bounds."; + } + + // Iterate over the elements in the row that are not on the diagonal and set them to zero. + T* valuePtr = valueStorage.data() + rowIndications[row]; + T* valuePtrEnd = valueStorage.data() + rowIndications[row + 1]; + uint_fast64_t* columnPtr = columnIndications.data() + rowIndications[row]; + + // 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 (valuePtr >= valuePtrEnd) { + throw storm::exceptions::InvalidStateException() << "Illegal call to SparseMatrix::makeRowAbsorbing: cannot make row " << row << " absorbing, but there is no entry in this row."; + } + + // If there is at least one entry in this row, we can just set it to one, modify its column value to the + // one given by the parameter and set all subsequent elements of this row to zero. + *valuePtr = storm::utility::constantOne(); + *columnPtr = column; + for (; valuePtr != valuePtrEnd; ++valuePtr) { + *valuePtr = storm::utility::constantZero(); + *columnPtr = 0; + } + } + + 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]); + } + } - for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { - if (columnConstraint.get(columnIndications[j])) { - ++subNonZeroEntries; - - if (index == columnIndications[j]) { - foundDiagonalElement = true; + ++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, bool insertDiagonalEntries) const { + return getSubmatrix(rowGroupConstraint, rowGroupConstraint, rowGroupIndices, insertDiagonalEntries); + } + + template + SparseMatrix SparseMatrix::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector const& rowGroupIndices, bool insertDiagonalEntries) 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) { + bool foundDiagonalElement = false; + + for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { + if (columnConstraint.get(columnIndications[j])) { + ++subNonZeroEntries; + + if (index == columnIndications[j]) { + foundDiagonalElement = true; + } } - } - } + } + + if (insertDiagonalEntries && !foundDiagonalElement) { + ++subNonZeroEntries; + } + } + } + + LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << rowGroupConstraint.getNumberOfSetBits() << "."); + + // Create and initialize resulting matrix. + SparseMatrix result(subRowCount, columnConstraint.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; + + // If we are requested to add missing diagonal entries, we need to make sure the corresponding rows + storm::storage::BitVector columnBitCountConstraint = columnConstraint; + if (insertDiagonalEntries) { + columnBitCountConstraint |= rowGroupConstraint; + } + for (auto index : columnBitCountConstraint) { + 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) { + bool insertedDiagonalElement = false; + + for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { + if (columnConstraint.get(columnIndications[j])) { + if (index == columnIndications[j]) { + insertedDiagonalElement = true; + } else if (insertDiagonalEntries && !insertedDiagonalElement && columnIndications[j] > index) { + result.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::constantZero()); + insertedDiagonalElement = true; + } + result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]); + } + } + if (insertDiagonalEntries && !insertedDiagonalElement) { + result.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::constantZero()); + } + + ++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 " << subRowCount << "x" << rowGroupConstraint.getNumberOfSetBits() << "."); - - // Create and initialize resulting matrix. - SparseMatrix result(subRowCount, columnConstraint.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; - - // If we are requested to add missing diagonal entries, we need to make sure the corresponding rows - storm::storage::BitVector columnBitCountConstraint = columnConstraint; - if (insertDiagonalEntries) { - columnBitCountConstraint |= rowGroupConstraint; - } - for (auto index : columnBitCountConstraint) { - 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) { - bool insertedDiagonalElement = false; + } + + 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]; - for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { - if (columnConstraint.get(columnIndications[j])) { - if (index == columnIndications[j]) { - insertedDiagonalElement = true; - } else if (insertDiagonalEntries && !insertedDiagonalElement && columnIndications[j] > index) { - result.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::constantZero()); - insertedDiagonalElement = true; - } - result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]); - } - } + // 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::constantZero()); + insertedDiagonalElement = true; + } + submatrix.addNextValue(rowGroupIndex, columnIndications[i], valueStorage[i]); + } if (insertDiagonalEntries && !insertedDiagonalElement) { - result.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::constantZero()); + submatrix.addNextValue(rowGroupIndex, rowGroupIndex, storm::utility::constantZero()); } - - ++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; + } + + // Finalize created matrix and return result. + submatrix.finalize(); + LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); + return submatrix; + } + + template + void SparseMatrix::convertToEquationSystem() { + invertDiagonal(); + negateAllNonDiagonalElements(); + } + + template + SparseMatrix SparseMatrix::transpose() const { + + uint_fast64_t rowCount = this->colCount; + uint_fast64_t colCount = this->rowCount; + uint_fast64_t nonZeroEntryCount = this->nonZeroEntryCount; + + std::vector rowIndications(rowCount + 1); + std::vector columnIndications(nonZeroEntryCount); + std::vector values(nonZeroEntryCount, T()); + + // First, we need to count how many entries each column has. + for (uint_fast64_t i = 0; i < this->rowCount; ++i) { + typename storm::storage::SparseMatrix::Rows rows = this->getRow(i); + for (auto const& transition : rows) { + if (transition.value() > 0) { + ++rowIndications[transition.column() + 1]; + } } - ++subNonZeroEntries; } - if (insertDiagonalEntries && !foundDiagonalElement) { - ++subNonZeroEntries; + + // Now compute the accumulated offsets. + for (uint_fast64_t i = 1; i < rowCount + 1; ++i) { + rowIndications[i] = rowIndications[i - 1] + rowIndications[i]; + } + + // Create an array that stores the index for the next value to be added for + // each row in the transposed matrix. Initially this corresponds to the previously + // computed accumulated offsets. + std::vector nextIndices = rowIndications; + + // Now we are ready to actually fill in the values of the transposed matrix. + for (uint_fast64_t i = 0; i < this->rowCount; ++i) { + typename storm::storage::SparseMatrix::Rows rows = this->getRow(i); + for (auto& transition : rows) { + if (transition.value() > 0) { + values[nextIndices[transition.column()]] = transition.value(); + columnIndications[nextIndices[transition.column()]++] = i; + } + } } + + storm::storage::SparseMatrix transposedMatrix(rowCount, colCount, + nonZeroEntryCount, + std::move(rowIndications), + std::move(columnIndications), + std::move(values)); + + return transposedMatrix; } - - 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::constantZero()); - insertedDiagonalElement = true; + + 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::constantOne(); + 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; } - submatrix.addNextValue(rowGroupIndex, columnIndications[i], valueStorage[i]); + + // 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!"; } - if (insertDiagonalEntries && !insertedDiagonalElement) { - submatrix.addNextValue(rowGroupIndex, rowGroupIndex, storm::utility::constantZero()); + + // 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; + } } } - - // Finalize created matrix and return result. - submatrix.finalize(); - LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); - return submatrix; - } - - template - void SparseMatrix::convertToEquationSystem() { - invertDiagonal(); - negateAllNonDiagonalElements(); - } - - template - SparseMatrix SparseMatrix::transpose() const { - - uint_fast64_t rowCount = this->colCount; - uint_fast64_t colCount = this->rowCount; - uint_fast64_t nonZeroEntryCount = this->nonZeroEntryCount; - - std::vector rowIndications(rowCount + 1); - std::vector columnIndications(nonZeroEntryCount); - std::vector values(nonZeroEntryCount, T()); - - // First, we need to count how many entries each column has. - for (uint_fast64_t i = 0; i < this->rowCount; ++i) { - typename storm::storage::SparseMatrix::Rows rows = this->getRow(i); - for (auto const& transition : rows) { - if (transition.value() > 0) { - ++rowIndications[transition.column() + 1]; - } - } - } - - // Now compute the accumulated offsets. - for (uint_fast64_t i = 1; i < rowCount + 1; ++i) { - rowIndications[i] = rowIndications[i - 1] + rowIndications[i]; - } - - // Create an array that stores the index for the next value to be added for - // each row in the transposed matrix. Initially this corresponds to the previously - // computed accumulated offsets. - std::vector nextIndices = rowIndications; - - // Now we are ready to actually fill in the values of the transposed matrix. - for (uint_fast64_t i = 0; i < this->rowCount; ++i) { - typename storm::storage::SparseMatrix::Rows rows = this->getRow(i); - for (auto& transition : rows) { - if (transition.value() > 0) { - values[nextIndices[transition.column()]] = transition.value(); - columnIndications[nextIndices[transition.column()]++] = i; - } - } - } - - storm::storage::SparseMatrix transposedMatrix(rowCount, colCount, - nonZeroEntryCount, - std::move(rowIndications), - std::move(columnIndications), - std::move(values)); - - return transposedMatrix; - } - - 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::constantOne(); - 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 std::pair, storm::storage::SparseMatrix> 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::constantOne(); - - // 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::constantZero(); - } - 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::constantZero()); - - // 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 { + + template + typename std::pair, storm::storage::SparseMatrix> 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::constantOne(); + + // 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::constantZero(); + } + 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::constantZero()); + + // 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)); + 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::constantZero(); - - for (auto elementIt = rowIt.begin(), elementIte = rowIt.end(); elementIt != elementIte; ++elementIt) { - *resultIt += elementIt.value() * vector[elementIt.column()]; + ConstRowIterator rowIt = this->begin(); + + for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++rowIt) { + *resultIt = storm::utility::constantZero(); + + 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::MutableRows SparseMatrix::getMutableRows(uint_fast64_t startRow, uint_fast64_t endRow) { - return MutableRows(this->valueStorage.data() + this->rowIndications[startRow], this->columnIndications.data() + this->rowIndications[startRow], this->rowIndications[endRow + 1] - this->rowIndications[startRow]); - } - - template - typename SparseMatrix::MutableRows SparseMatrix::getMutableRow(uint_fast64_t row) { - return getMutableRows(row, row); - } - - 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 - typename SparseMatrix::ValueIterator SparseMatrix::valueIteratorBegin(uint_fast64_t row) { - return &(this->valueStorage[0]) + this->rowIndications[row]; - } - - template - typename SparseMatrix::ValueIterator SparseMatrix::valueIteratorEnd(uint_fast64_t row) { - return &(this->valueStorage[0]) + this->rowIndications[row + 1]; - } - - template - T SparseMatrix::getRowSum(uint_fast64_t row) const { - T sum = storm::utility::constantZero(); - 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; + } + + 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::MutableRows SparseMatrix::getMutableRows(uint_fast64_t startRow, uint_fast64_t endRow) { + return MutableRows(this->valueStorage.data() + this->rowIndications[startRow], this->columnIndications.data() + this->rowIndications[startRow], this->rowIndications[endRow + 1] - this->rowIndications[startRow]); + } + + template + typename SparseMatrix::MutableRows SparseMatrix::getMutableRow(uint_fast64_t row) { + return getMutableRows(row, row); + } + + 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 + typename SparseMatrix::ValueIterator SparseMatrix::valueIteratorBegin(uint_fast64_t row) { + return &(this->valueStorage[0]) + this->rowIndications[row]; + } + + template + typename SparseMatrix::ValueIterator SparseMatrix::valueIteratorEnd(uint_fast64_t row) { + return &(this->valueStorage[0]) + this->rowIndications[row + 1]; + } + + template + T SparseMatrix::getRowSum(uint_fast64_t row) const { + T sum = storm::utility::constantZero(); + 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::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 << std::setprecision(8) << 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, boost::hash_range(valueStorage.begin(), valueStorage.end())); - boost::hash_combine(result, boost::hash_range(columnIndications.begin(), columnIndications.end())); - boost::hash_combine(result, boost::hash_range(rowIndications.begin(), rowIndications.end())); - - 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::constantZero()); - 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. - + } + return true; + } + + 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 << std::setprecision(8) << 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, boost::hash_range(valueStorage.begin(), valueStorage.end())); + boost::hash_combine(result, boost::hash_range(columnIndications.begin(), columnIndications.end())); + boost::hash_combine(result, boost::hash_range(rowIndications.begin(), rowIndications.end())); + + 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::constantZero()); + 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::constantZero(); - - 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>; - + + 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::constantZero(); + + 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 storage } // namespace storm diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 7970ae3f2..ab1e0d097 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -154,8 +154,7 @@ namespace storm { /*! * This class represents a number of consecutive rows of the matrix. */ - template - class BaseRows { + class rows { public: /*! * Constructs an object that represents the rows defined by the value of the first entry, the column @@ -165,25 +164,65 @@ namespace storm { * @param columnPtr A pointer to the column of the first entry of the rows. * @param entryCount The number of entrys in the rows. */ - Rows(ValueType* valuePtr, uint_fast64_t const* columnPtr, uint_fast64_t entryCount); + rows(T* 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. */ - BaseIterator begin(); + iterator begin(); /*! * Retrieves an iterator that points past the last entry of the rows. * * @return An iterator that points past the last entry of the rows. */ - BaseIterator end() const; + iterator end(); private: // The pointer to the value of the first entry. - ValueType* valuePtr; + T* valuePtr; + + // The pointer to the column of the first entry. + uint_fast64_t const* columnPtr; + + // The number of non-zero entries in the rows. + uint_fast64_t entryCount; + }; + + /*! + * This class represents a number of consecutive rows of the matrix. + */ + class const_rows { + public: + /*! + * Constructs an object that represents the rows defined by the value of the first entry, the column + * of the first entry and the number of entries in this row set. + * + * @param valuePtr A pointer to the value of the first entry of the rows. + * @param columnPtr A pointer to the column of the first entry of the rows. + * @param entryCount The number of entrys in the rows. + */ + const_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. + */ + const_iterator begin() const; + + /*! + * Retrieves an iterator that points past the last entry of the rows. + * + * @return An iterator that points past the last entry of the rows. + */ + const_iterator end() const; + + private: + // The pointer to the value of the first entry. + T const* valuePtr; // The pointer to the column of the first entry. uint_fast64_t const* columnPtr; @@ -191,9 +230,6 @@ namespace storm { // The number of non-zero entries in the rows. uint_fast64_t entryCount; }; - - typedef BaseRows rows; - typedef BaseRows const_rows; /*! * An enum representing the internal state of the matrix. After creation, the matrix is UNINITIALIZED. @@ -217,7 +253,7 @@ namespace storm { * @param size The number of rows and columns of the matrix. * @param entries The number of entries of the matrix. */ - SparseMatrix(uint_fast64_t size = 0, uint_fast64_t entries); + SparseMatrix(uint_fast64_t size = 0, uint_fast64_t entries = 0); /*! * Constructs a sparse matrix by performing a deep-copy of the given matrix. @@ -236,27 +272,26 @@ namespace storm { /*! * Constructs a sparse matrix by moving the given contents. * - * @param colCount The number of columns of the matrix. + * @param columnCount The number of columns of the matrix. * @param rowIndications The row indications vector of the matrix to be constructed. * @param columnIndications The column indications vector of the matrix to be constructed. * @param values The vector containing the values of the entries in the matrix. */ - SparseMatrix(uint_fast64_t colCount, std::vector&& rowIndications, - std::vector&& columnIndications, std::vector&& values); + SparseMatrix(uint_fast64_t columnCount, std::vector&& rowIndications, std::vector&& columnIndications, std::vector&& values); /*! * Assigns the contents of the given matrix to the current one by deep-copying its contents. * * @param other The matrix from which to copy-assign. */ - storm::storage::SparseMatrix& operator=(SparseMatrix const& other); + SparseMatrix& operator=(SparseMatrix const& other); /*! * Assigns the contents of the given matrix to the current one by moving its contents. * * @param other The matrix from which to move to contents. */ - storm::storage::SparseMatrix& operator=(SparseMatrix&& other); + SparseMatrix& operator=(SparseMatrix&& other); /*! * Sets the matrix entry at the given row and column to the given value. After all entries have been added, @@ -268,15 +303,10 @@ namespace storm { * treated as empty. If these constraints are not met, an exception is thrown. * * @param row The row in which the matrix entry is to be set. - * @param col The column in which the matrix entry is to be set. + * @param column The column in which the matrix entry is to be set. * @param value The value that is to be set at the specified row and column. */ - void addNextValue(uint_fast64_t row, uint_fast64_t col, T const& value); - - /*! - * Inserts an empty row in the matrix. - */ - void insertEmptyRow(); + void addNextValue(uint_fast64_t row, uint_fast64_t column, T const& value); /* * Finalizes the sparse matrix to indicate that initialization process has been completed and the matrix @@ -631,8 +661,8 @@ namespace storm { #ifdef STORM_HAVE_INTELTBB /*! - * This function is a helper for Parallel Execution of the multipliyWithVector functionality. - * It uses Intels TBB parallel_for paradigm to split up the row/vector multiplication and summation + * This function is a helper for the parallel execution of the multipliyWithVector method. + * It uses Intel's TBB parallel_for paradigm to split up the row/vector multiplication and summation. */ template class tbbHelper_MatrixRowVectorScalarProduct { @@ -645,7 +675,6 @@ namespace storm { V * resultVector; V const* vectorX; M const* matrixA; - }; #endif