sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Further work on sparse matrix implementation. 

Former-commit-id: df4eb9c476
tempestpy_adaptions
dehnert 11 years ago
parent
4376708a46
commit
ef041982b5
2 changed files with 56 additions and 74 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 111
      src/storage/SparseMatrix.cpp
  2. 19
      src/storage/SparseMatrix.h

111
src/storage/SparseMatrix.cpp
View File

@ -323,12 +323,12 @@ namespace storm {
} }
template<typename T> template<typename T>
std::vector<T> SparseMatrix<T>::getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, storm::storage::BitVector const& columnConstraint, uint_fast64_t numberOfRows) const {
std::vector<T> result(numberOfRows);
uint_fast64_t currentRowCount = 0;
std::vector<T> SparseMatrix<T>::getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, storm::storage::BitVector const& columnConstraint) const {
std::vector<T> result;
result.reserve(rowGroupConstraint.getNumberOfSetBits());
for (auto rowGroup : rowGroupConstraint) { for (auto rowGroup : rowGroupConstraint) {
for (uint_fast64_t row = rowGroupIndices[rowGroup]; row < rowGroupIndices[rowGroup + 1]; ++row) { 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; return result;
@ -336,33 +336,28 @@ namespace storm {
template<typename T> template<typename T>
SparseMatrix<T> SparseMatrix<T>::getSubmatrix(storm::storage::BitVector const& constraint) const { SparseMatrix<T> SparseMatrix<T>::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) { 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 (auto rowIndex : constraint) {
for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) { for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) {
if (constraint.get(columnIndications[i])) { if (constraint.get(columnIndications[i])) {
++subNonZeroEntries;
++subEntries;
} }
} }
} }
// Create and initialize resulting matrix. // 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<uint_fast64_t> bitsSetBeforeIndex; std::vector<uint_fast64_t> bitsSetBeforeIndex;
bitsSetBeforeIndex.reserve(colCount);
bitsSetBeforeIndex.reserve(columnCount);
// Compute the information to fill this vector. // Compute the information to fill this vector.
uint_fast64_t lastIndex = 0; uint_fast64_t lastIndex = 0;
@ -375,7 +370,7 @@ namespace storm {
++currentNumberOfSetBits; ++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; uint_fast64_t rowCount = 0;
for (auto rowIndex : constraint) { for (auto rowIndex : constraint) {
for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) { for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) {
@ -387,9 +382,8 @@ namespace storm {
++rowCount; ++rowCount;
} }
// Finalize sub-matrix and return result.
// Finalize submatrix and return result.
result.finalize(); result.finalize();
LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix.");
return result; return result;
} }
@ -400,19 +394,17 @@ namespace storm {
template<typename T> template<typename T>
SparseMatrix<T> SparseMatrix<T>::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries) const { SparseMatrix<T> SparseMatrix<T>::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector<uint_fast64_t> 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) { 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) { for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) {
bool foundDiagonalElement = false; bool foundDiagonalElement = false;
for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) {
if (columnConstraint.get(columnIndications[j])) { if (columnConstraint.get(columnIndications[j])) {
++subNonZeroEntries;
++subEntries;
if (index == columnIndications[j]) { if (index == columnIndications[j]) {
foundDiagonalElement = true; 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) { if (insertDiagonalEntries && !foundDiagonalElement) {
++subNonZeroEntries;
++subEntries;
} }
} }
} }
LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << rowGroupConstraint.getNumberOfSetBits() << ".");
// Create and initialize resulting matrix. // 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<uint_fast64_t> bitsSetBeforeIndex; std::vector<uint_fast64_t> bitsSetBeforeIndex;
bitsSetBeforeIndex.reserve(colCount);
bitsSetBeforeIndex.reserve(columnCount);
// Compute the information to fill this vector. // Compute the information to fill this vector.
uint_fast64_t lastIndex = 0; uint_fast64_t lastIndex = 0;
uint_fast64_t currentNumberOfSetBits = 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; storm::storage::BitVector columnBitCountConstraint = columnConstraint;
if (insertDiagonalEntries) { if (insertDiagonalEntries) {
columnBitCountConstraint |= rowGroupConstraint; columnBitCountConstraint |= rowGroupConstraint;
@ -479,19 +470,15 @@ namespace storm {
} }
} }
// Finalize sub-matrix and return result.
result.finalize(); result.finalize();
LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix.");
return result; return result;
} }
template<typename T> template<typename T>
SparseMatrix<T> SparseMatrix<T>::getSubmatrix(std::vector<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries) const { SparseMatrix<T> SparseMatrix<T>::getSubmatrix(std::vector<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> 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) { for (uint_fast64_t rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) {
// Determine which row we need to select from the current row group. // Determine which row we need to select from the current row group.
uint_fast64_t rowToCopy = rowGroupIndices[rowGroupIndex] + rowGroupToRowIndexMapping[rowGroupIndex]; uint_fast64_t rowToCopy = rowGroupIndices[rowGroupIndex] + rowGroupToRowIndexMapping[rowGroupIndex];
@ -502,26 +489,23 @@ namespace storm {
if (columnIndications[i] == rowGroupIndex) { if (columnIndications[i] == rowGroupIndex) {
foundDiagonalElement = true; foundDiagonalElement = true;
} }
++subNonZeroEntries;
++subEntries;
} }
if (insertDiagonalEntries && !foundDiagonalElement) { 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. // Now create the matrix to be returned with the appropriate size.
SparseMatrix<T> submatrix(rowGroupIndices.size() - 1, colCount);
submatrix.initialize(subNonZeroEntries);
SparseMatrix<T> submatrix(rowGroupIndices.size() - 1, columnCount, subEntries);
// Copy over the selected lines from the source matrix. // Copy over the selected lines from the source matrix.
for (uint_fast64_t rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) { for (uint_fast64_t rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) {
// Determine which row we need to select from the current row group. // Determine which row we need to select from the current row group.
uint_fast64_t rowToCopy = rowGroupIndices[rowGroupIndex] + rowGroupToRowIndexMapping[rowGroupIndex]; 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; bool insertedDiagonalElement = false;
for (uint_fast64_t i = rowIndications[rowToCopy], rowEnd = rowIndications[rowToCopy + 1]; i < rowEnd; ++i) { for (uint_fast64_t i = rowIndications[rowToCopy], rowEnd = rowIndications[rowToCopy + 1]; i < rowEnd; ++i) {
if (columnIndications[i] == rowGroupIndex) { if (columnIndications[i] == rowGroupIndex) {
@ -539,26 +523,19 @@ namespace storm {
// Finalize created matrix and return result. // Finalize created matrix and return result.
submatrix.finalize(); submatrix.finalize();
LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix.");
return submatrix; return submatrix;
} }
template<typename T>
void SparseMatrix<T>::convertToEquationSystem() {
invertDiagonal();
negateAllNonDiagonalElements();
}
template <typename T> template <typename T>
SparseMatrix<T> SparseMatrix<T>::transpose() const { SparseMatrix<T> SparseMatrix<T>::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<uint_fast64_t> rowIndications(rowCount + 1); std::vector<uint_fast64_t> rowIndications(rowCount + 1);
std::vector<uint_fast64_t> columnIndications(nonZeroEntryCount);
std::vector<T> values(nonZeroEntryCount, T());
std::vector<uint_fast64_t> columnIndications(entryCount);
std::vector<T> values(entryCount);
// First, we need to count how many entries each column has. // First, we need to count how many entries each column has.
for (uint_fast64_t i = 0; i < this->rowCount; ++i) { for (uint_fast64_t i = 0; i < this->rowCount; ++i) {
@ -600,6 +577,12 @@ namespace storm {
return transposedMatrix; return transposedMatrix;
} }
template<typename T>
void SparseMatrix<T>::convertToEquationSystem() {
invertDiagonal();
negateAllNonDiagonalEntries();
}
template<typename T> template<typename T>
void SparseMatrix<T>::invertDiagonal() { void SparseMatrix<T>::invertDiagonal() {
// Check if the matrix is square, because only then it makes sense to perform this // Check if the matrix is square, because only then it makes sense to perform this

19
src/storage/SparseMatrix.h
View File

@ -376,9 +376,8 @@ namespace storm {
T getConstrainedRowSum(uint_fast64_t row, storm::storage::BitVector const& columns) const; 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 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. * @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 * @return A vector whose entries represent the sums of selected columns for all rows in selected row
* groups. * groups.
*/ */
std::vector<T> getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices, storm::storage::BitVector const& columnConstraint, uint_fast64_t numberOfRows) const;
std::vector<T> getConstrainedRowSumVector(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> 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 * 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<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> const& rowGroupIndices, bool insertDiagonalEntries = true) const; SparseMatrix getSubmatrix(std::vector<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> 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. * Transposes the matrix.
* *
@ -458,6 +452,11 @@ namespace storm {
*/ */
storm::storage::SparseMatrix<T> transpose() const; storm::storage::SparseMatrix<T> 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. * 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. * 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. * 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. * Calculates the Jacobi decomposition of this sparse matrix. For this operation, the matrix must be square.

Write Preview
Loading…
Cancel
Save