diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index f2ffb95b3..97b0849ea 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -323,12 +323,12 @@ namespace storm { } 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; + std::vector SparseMatrix::getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices, storm::storage::BitVector const& columnConstraint) const { + std::vector result; + result.reserve(rowGroupConstraint.getNumberOfSetBits()); for (auto rowGroup : rowGroupConstraint) { for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) { - result[currentRowCount++] = getConstrainedRowSum(row, columnConstraint); + result.push_back(getConstrainedRowSum(row, columnConstraint)); } } return result; @@ -336,33 +336,28 @@ namespace storm { 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. + // Check whether we select at least some rows and columns. 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."); + throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrix::getSubmatrix: cannot create empty submatrix."; } - // First, we need to determine the number of non-zero entries of the - // sub-matrix. - uint_fast64_t subNonZeroEntries = 0; + // First, we need to determine the number of entries of the submatrix. + uint_fast64_t subEntries = 0; for (auto rowIndex : constraint) { for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) { if (constraint.get(columnIndications[i])) { - ++subNonZeroEntries; + ++subEntries; } } } // Create and initialize resulting matrix. - SparseMatrix result(constraint.getNumberOfSetBits()); - result.initialize(subNonZeroEntries); + SparseMatrix result(constraint.getNumberOfSetBits(), subEntries); - // 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. + // 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); + bitsSetBeforeIndex.reserve(columnCount); // Compute the information to fill this vector. uint_fast64_t lastIndex = 0; @@ -375,7 +370,7 @@ namespace storm { ++currentNumberOfSetBits; } - // Copy over selected entries. + // Copy over selected entries and use the previously computed vector to get the column offset. uint_fast64_t rowCount = 0; for (auto rowIndex : constraint) { for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) { @@ -387,9 +382,8 @@ namespace storm { ++rowCount; } - // Finalize sub-matrix and return result. + // Finalize submatrix and return result. result.finalize(); - LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); return result; } @@ -400,19 +394,17 @@ namespace storm { 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; + // First, we need to determine the number of entries and the number of rows of the submatrix. + uint_fast64_t subEntries = 0; + uint_fast64_t subRows = 0; for (auto index : rowGroupConstraint) { - subRowCount += rowGroupIndices[index + 1] - rowGroupIndices[index]; + subRows += 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; + ++subEntries; if (index == columnIndications[j]) { foundDiagonalElement = true; @@ -420,28 +412,27 @@ namespace storm { } } + // If requested, we need to reserve one entry more for inserting the diagonal zero entry. if (insertDiagonalEntries && !foundDiagonalElement) { - ++subNonZeroEntries; + ++subEntries; } } } - 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); + SparseMatrix result(subRows, columnConstraint.getNumberOfSetBits(), subEntries); - // 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. + // 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); + bitsSetBeforeIndex.reserve(columnCount); // 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 + // If we are requested to add missing diagonal entries, we need to make sure the corresponding rows are also + // taken. storm::storage::BitVector columnBitCountConstraint = columnConstraint; if (insertDiagonalEntries) { columnBitCountConstraint |= rowGroupConstraint; @@ -479,19 +470,15 @@ namespace storm { } } - // 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; + // First, we need to count how many non-zero entries the resulting matrix will have and reserve space for + // diagonal entries if requested. + uint_fast64_t subEntries = 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]; @@ -502,26 +489,23 @@ namespace storm { if (columnIndications[i] == rowGroupIndex) { foundDiagonalElement = true; } - ++subNonZeroEntries; + ++subEntries; } if (insertDiagonalEntries && !foundDiagonalElement) { - ++subNonZeroEntries; + ++subEntries; } } - 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); + SparseMatrix submatrix(rowGroupIndices.size() - 1, columnCount, subEntries); // 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. + // 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) { @@ -539,26 +523,19 @@ namespace storm { // 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; + uint_fast64_t rowCount = this->columnCount; + uint_fast64_t columnCount = this->rowCount; + uint_fast64_t entryCount = this->nonZeroEntryCount; std::vector rowIndications(rowCount + 1); - std::vector columnIndications(nonZeroEntryCount); - std::vector values(nonZeroEntryCount, T()); + std::vector columnIndications(entryCount); + std::vector values(entryCount); // First, we need to count how many entries each column has. for (uint_fast64_t i = 0; i < this->rowCount; ++i) { @@ -600,6 +577,12 @@ namespace storm { return transposedMatrix; } + template + void SparseMatrix::convertToEquationSystem() { + invertDiagonal(); + negateAllNonDiagonalEntries(); + } + template void SparseMatrix::invertDiagonal() { // Check if the matrix is square, because only then it makes sense to perform this diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index ab1e0d097..36413c08a 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -376,9 +376,8 @@ namespace storm { T getConstrainedRowSum(uint_fast64_t row, storm::storage::BitVector const& columns) const; /*! - * Computes a vector whose i-th entry is the sum of the entries in the i-th row where only those entries are - * added that are in selected columns. Likewise, this sum is only computed for selected rows and set to zero - * for all unselected rows. + * Computes a vector whose i-th entry is the sum of the entries in the i-th selected row where only those + * entries are added that are in selected columns. * * @param rowConstraint A bit vector that indicates for which rows to compute the constrained sum. * @param columnConstraint A bit vector that indicates which columns to add in the selected rows. @@ -397,7 +396,7 @@ namespace storm { * @return A vector whose entries represent the sums of selected columns for all rows in selected row * groups. */ - std::vector getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices, storm::storage::BitVector const& columnConstraint, uint_fast64_t numberOfRows) const; + std::vector getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices, storm::storage::BitVector const& columnConstraint) const; /*! * Creates a submatrix of the current matrix by dropping all rows and columns whose bits are not @@ -446,11 +445,6 @@ namespace storm { */ SparseMatrix getSubmatrix(std::vector const& rowGroupToRowIndexMapping, std::vector const& rowGroupIndices, bool insertDiagonalEntries = true) const; - /*! - * Transforms the matrix into an equation system. That is, it transforms the matrix A into a matrix (1-A). - */ - void convertToEquationSystem(); - /*! * Transposes the matrix. * @@ -458,6 +452,11 @@ namespace storm { */ storm::storage::SparseMatrix transpose() const; + /*! + * Transforms the matrix into an equation system. That is, it transforms the matrix A into a matrix (1-A). + */ + void convertToEquationSystem(); + /*! * Inverts all entries on the diagonal, i.e. sets the diagonal values to one minus their previous value. * Requires the matrix to contain each diagonal entry and to be square. @@ -467,7 +466,7 @@ namespace storm { /*! * Negates (w.r.t. addition) all entries that are not on the diagonal. */ - void negateAllNonDiagonalentries(); + void negateAllNonDiagonalEntries(); /*! * Calculates the Jacobi decomposition of this sparse matrix. For this operation, the matrix must be square.