From dd035f7f5eb23ab492d8777ef93e347348a36a8c Mon Sep 17 00:00:00 2001 From: dehnert Date: Fri, 8 Sep 2017 13:57:57 +0200 Subject: [PATCH] allow for summand in matrix-vector multiplication --- .../IterativeMinMaxLinearEquationSolver.cpp | 1 - .../solver/NativeLinearEquationSolver.cpp | 15 +--- src/storm/storage/SparseMatrix.cpp | 74 +++++++++++-------- src/storm/storage/SparseMatrix.h | 10 ++- .../NativeMinMaxLinearEquationSolverTest.cpp | 1 + 5 files changed, 56 insertions(+), 45 deletions(-) diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 314c97c2b..f27bd97bf 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -256,7 +256,6 @@ namespace storm { storm::utility::vector::selectVectorValues(*auxiliaryRowGroupVector, this->getInitialScheduler(), this->A->getRowGroupIndices(), b); // Solve the resulting equation system. - // Note that the linEqSolver might consider a slightly different interpretation of "equalModuloPrecision". Hence, we iteratively increase its precision. auto submatrixSolver = this->linearEquationSolverFactory->create(std::move(submatrix)); submatrixSolver->setCachingEnabled(true); if (this->lowerBound) { submatrixSolver->setLowerBound(this->lowerBound.get()); } diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 339e56c6f..ae4824672 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -207,22 +207,15 @@ namespace storm { template void NativeLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { if (&x != &result) { - A->multiplyWithVector(x, result); - if (b != nullptr) { - storm::utility::vector::addVectors(result, *b, result); - } + A->multiplyWithVector(x, result, b); } else { // If the two vectors are aliases, we need to create a temporary. - if(!this->cachedRowVector) { + if (!this->cachedRowVector) { this->cachedRowVector = std::make_unique>(getMatrixRowCount()); } - A->multiplyWithVector(x, *this->cachedRowVector); - if (b != nullptr) { - storm::utility::vector::addVectors(*this->cachedRowVector, *b, result); - } else { - result.swap(*this->cachedRowVector); - } + A->multiplyWithVector(x, *this->cachedRowVector, b); + result.swap(*this->cachedRowVector); if (!this->isCachingEnabled()) { clearCache(); diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index 49d0cddad..2d422e0c7 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -76,7 +76,7 @@ namespace storm { template bool MatrixEntry::operator!=(MatrixEntry const& other) const { - return !(*this == other); + return !(*this == other); } template @@ -210,7 +210,7 @@ namespace storm { bool hasEntries = currentEntryCount != 0; uint_fast64_t rowCount = hasEntries ? lastRow + 1 : 0; - + // If the last row group was empty, we need to add one more to the row count, because otherwise this empty row is not counted. if (hasCustomRowGrouping) { if (lastRow < rowGroupIndices->back()) { @@ -304,7 +304,7 @@ namespace storm { template void SparseMatrixBuilder::replaceColumns(std::vector const& replacements, index_type offset) { index_type maxColumn = 0; - + for (index_type row = 0; row < rowIndications.size(); ++row) { bool changed = false; auto startRow = std::next(columnsAndValues.begin(), rowIndications[row]); @@ -330,11 +330,11 @@ namespace storm { }), "Columns not sorted."); } } - + highestColumn = maxColumn; lastColumn = columnsAndValues.empty() ? 0 : columnsAndValues[columnsAndValues.size() - 1].getColumn(); } - + template SparseMatrix::rows::rows(iterator begin, index_type entryCount) : beginIterator(begin), entryCount(entryCount) { // Intentionally left empty. @@ -424,7 +424,6 @@ namespace storm { rowGroupIndices = other.rowGroupIndices; trivialRowGrouping = other.trivialRowGrouping; } - return *this; } @@ -442,7 +441,6 @@ namespace storm { rowGroupIndices = std::move(other.rowGroupIndices); trivialRowGrouping = other.trivialRowGrouping; } - return *this; } @@ -583,7 +581,7 @@ namespace storm { } return rowGroupIndices.get(); } - + template storm::storage::BitVector SparseMatrix::getRowFilter(storm::storage::BitVector const& groupConstraint) const { storm::storage::BitVector res(this->getRowCount(), false); @@ -1125,7 +1123,7 @@ namespace storm { return transposedMatrix; } - + template SparseMatrix SparseMatrix::transposeSelectedRowsFromRowGroups(std::vector const& rowGroupChoices, bool keepZeros) const { index_type rowCount = this->getColumnCount(); @@ -1291,22 +1289,22 @@ namespace storm { return result; } - + template - void SparseMatrix::multiplyWithVector(std::vector const& vector, std::vector& result) const { + void SparseMatrix::multiplyWithVector(std::vector const& vector, std::vector& result, std::vector const* summand) const { #ifdef STORM_HAVE_INTELTBB if (this->getNonzeroEntryCount() > 10000) { - return this->multiplyWithVectorParallel(vector, result); + return this->multiplyWithVectorParallel(vector, result, summand); } else { - return this->multiplyWithVectorSequential(vector, result); + return this->multiplyWithVectorSequential(vector, result, summand); } #else - return multiplyWithVectorSequential(vector, result); + return multiplyWithVectorSequential(vector, result, summand); #endif } - + template - void SparseMatrix::multiplyWithVectorSequential(std::vector const& vector, std::vector& result) const { + void SparseMatrix::multiplyWithVectorSequential(std::vector const& vector, std::vector& result, std::vector const* summand) const { if (&vector == &result) { STORM_LOG_WARN("Matrix-vector-multiplication invoked but the target vector uses the same memory as the input vector. This requires to allocate auxiliary memory."); std::vector tmpVector(this->getRowCount()); @@ -1318,20 +1316,29 @@ namespace storm { std::vector::const_iterator rowIterator = rowIndications.begin(); typename std::vector::iterator resultIterator = result.begin(); typename std::vector::iterator resultIteratorEnd = result.end(); - + typename std::vector::const_iterator summandIterator; + if (summand) { + summandIterator = summand->begin(); + } + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) { - *resultIterator = storm::utility::zero(); - + if (summand) { + *resultIterator = *summandIterator; + ++summandIterator; + } else { + *resultIterator = storm::utility::zero(); + } + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { *resultIterator += it->getValue() * vector[it->getColumn()]; } } } } - + #ifdef STORM_HAVE_INTELTBB template - void SparseMatrix::multiplyWithVectorParallel(std::vector const& vector, std::vector& result) const { + void SparseMatrix::multiplyWithVectorParallel(std::vector const& vector, std::vector& result, std::vector const* summand) const { if (&vector == &result) { STORM_LOG_WARN("Matrix-vector-multiplication invoked but the target vector uses the same memory as the input vector. This requires to allocate auxiliary memory."); std::vector tmpVector(this->getRowCount()); @@ -1348,10 +1355,19 @@ namespace storm { std::vector::const_iterator rowIteratorEnd = this->rowIndications.begin() + endRow; typename std::vector::iterator resultIterator = result.begin() + startRow; typename std::vector::iterator resultIteratorEnd = result.begin() + endRow; - + typename std::vector::const_iterator summandIterator; + if (summand) { + summandIterator = summand->begin() + startRow; + } + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) { - *resultIterator = storm::utility::zero(); - + if (summand) { + *resultIterator = *summandIterator; + ++summandIterator; + } else { + *resultIterator = storm::utility::zero(); + } + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { *resultIterator += it->getValue() * vector[it->getColumn()]; } @@ -1433,7 +1449,7 @@ namespace storm { ++row; } } - + template void SparseMatrix::divideRowsInPlace(std::vector const& divisors) { STORM_LOG_ASSERT(divisors.size() == this->getRowCount(), "Can not divide rows: Number of rows and number of divisors do not match."); @@ -1751,7 +1767,7 @@ namespace storm { template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; - + template class MatrixEntry; template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); @@ -1792,7 +1808,7 @@ namespace storm { template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; #endif - + #if defined(STORM_HAVE_GMP) template class MatrixEntry::index_type, GmpRationalNumber>; template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); @@ -1802,7 +1818,7 @@ namespace storm { template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; #endif - + // Rational Function template class MatrixEntry::index_type, RationalFunction>; template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); @@ -1814,7 +1830,7 @@ namespace storm { template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; - + // Intervals template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template class MatrixEntry::index_type, Interval>; diff --git a/src/storm/storage/SparseMatrix.h b/src/storm/storage/SparseMatrix.h index 6f4dba96f..50712ced6 100644 --- a/src/storm/storage/SparseMatrix.h +++ b/src/storm/storage/SparseMatrix.h @@ -773,9 +773,10 @@ namespace storm { * * @param vector The vector with which to multiply the matrix. * @param result The vector that is supposed to hold the result of the multiplication after the operation. + * @param summand If given, this summand will be added to the result of the multiplication. * @return The product of the matrix and the given vector as the content of the given result vector. */ - void multiplyWithVector(std::vector const& vector, std::vector& result) const; + void multiplyWithVector(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; /*! * Multiplies a single row of the matrix with the given vector and returns the result @@ -826,9 +827,10 @@ namespace storm { * * @param vector The vector with which to multiply the matrix. * @param result The vector that is supposed to hold the result of the multiplication after the operation. + * @param summand If given, this summand will be added to the result of the multiplication. * @return The product of the matrix and the given vector as the content of the given result vector. */ - void multiplyWithVectorSequential(std::vector const& vector, std::vector& result) const; + void multiplyWithVectorSequential(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; #ifdef STORM_HAVE_INTELTBB /*! @@ -837,9 +839,10 @@ namespace storm { * * @param vector The vector with which to multiply the matrix. * @param result The vector that is supposed to hold the result of the multiplication after the operation. + * @param summand If given, this summand will be added to the result. * @return The product of the matrix and the given vector as the content of the given result vector. */ - void multiplyWithVectorParallel(std::vector const& vector, std::vector& result) const; + void multiplyWithVectorParallel(std::vector const& vector, std::vector& result, std::vector const* summand = nullptr) const; #endif /*! @@ -1077,7 +1080,6 @@ namespace storm { // A vector indicating the row groups of the matrix. This needs to be mutible in case we create it on-the-fly. mutable boost::optional> rowGroupIndices; - }; #ifdef STORM_HAVE_CARL diff --git a/src/test/storm/solver/NativeMinMaxLinearEquationSolverTest.cpp b/src/test/storm/solver/NativeMinMaxLinearEquationSolverTest.cpp index de724e39f..bfd1defcc 100644 --- a/src/test/storm/solver/NativeMinMaxLinearEquationSolverTest.cpp +++ b/src/test/storm/solver/NativeMinMaxLinearEquationSolverTest.cpp @@ -22,6 +22,7 @@ TEST(NativeMinMaxLinearEquationSolver, SolveWithStandardOptions) { auto solver = factory.create(A); ASSERT_NO_THROW(solver->solveEquations(storm::OptimizationDirection::Minimize, x, b)); + ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::getModule().getPrecision()); ASSERT_NO_THROW(solver->solveEquations(storm::OptimizationDirection::Maximize, x, b));