diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 845e7495e..3b80a215b 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -1,4 +1,5 @@ #include +#include #include "storm/solver/IterativeMinMaxLinearEquationSolver.h" @@ -614,7 +615,12 @@ namespace storm { template class SoundValueIterationHelper { public: - SoundValueIterationHelper(std::vector& x, std::vector& y, bool relative, ValueType const& precision, uint64_t sizeOfLargestRowGroup) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), minIndex(0), maxIndex(0), relative(relative), precision(precision) { + + typedef uint32_t IndexType; + + SoundValueIterationHelper(storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector& y, bool relative, ValueType const& precision) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), minIndex(0), maxIndex(0), relative(relative), precision(precision), rowGroupIndices(matrix.getRowGroupIndices()) { + STORM_LOG_THROW(matrix.getEntryCount() < std::numeric_limits::max(), storm::exceptions::NotSupportedException, "The number of matrix entries is too large for the selected index type."); + uint64_t sizeOfLargestRowGroup = matrix.getSizeOfLargestRowGroup(); xTmp.resize(sizeOfLargestRowGroup); yTmp.resize(sizeOfLargestRowGroup); x.assign(x.size(), storm::utility::zero()); @@ -623,6 +629,22 @@ namespace storm { decisionValueBlocks = false; convergencePhase1 = true; firstIndexViolatingConvergence = 0; + + numRows = matrix.getRowCount(); + matrixValues.clear(); + matrixColumns.clear(); + rowIndications.clear(); + matrixValues.reserve(matrix.getNonzeroEntryCount()); + matrixColumns.reserve(matrix.getColumnCount()); + rowIndications.reserve(numRows + 1); + rowIndications.push_back(0); + for (IndexType r = 0; r < numRows; ++r) { + for (auto const& entry : matrix.getRow(r)) { + matrixValues.push_back(entry.getValue()); + matrixColumns.push_back(entry.getColumn()); + } + rowIndications.push_back(matrixValues.size()); + } } @@ -666,30 +688,38 @@ namespace storm { return maximize(dir) ? minIndex : maxIndex; } - void multiplyRow(uint64_t const& row, storm::storage::SparseMatrix const& A, storm::solver::Multiplier const& multiplier, ValueType const& bi, ValueType& xi, ValueType& yi) { + void multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi) { + assert(rowIndex < numRows); ValueType xRes = bi; ValueType yRes = storm::utility::zero(); - multiplier.multiplyRow2(row, x, xRes, y, yRes); + + auto entryIt = matrixValues.begin() + rowIndications[rowIndex]; + auto entryItE = matrixValues.begin() + rowIndications[rowIndex + 1]; + auto colIt = matrixColumns.begin() + rowIndications[rowIndex]; + for (; entryIt != entryItE; ++entryIt, ++colIt) { + xRes += *entryIt * x[*colIt]; + yRes += *entryIt * y[*colIt]; + } xi = std::move(xRes); yi = std::move(yRes); } template - void performIterationStep(storm::storage::SparseMatrix const& A, storm::solver::Multiplier const& multiplier, std::vector const& b) { + void performIterationStep(std::vector const& b) { if (!decisionValueBlocks) { - performIterationStepUpdateDecisionValue(A, multiplier, b); + performIterationStepUpdateDecisionValue(b); } else { assert(decisionValue == getPrimaryBound()); auto xIt = x.rbegin(); auto yIt = y.rbegin(); - auto groupStartIt = A.getRowGroupIndices().rbegin(); + auto groupStartIt = rowGroupIndices.rbegin(); uint64_t groupEnd = *groupStartIt; ++groupStartIt; - for (auto groupStartIte = A.getRowGroupIndices().rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { + for (auto groupStartIte = rowGroupIndices.rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { // Perform the iteration for the first row in the group - uint64_t row = *groupStartIt; + IndexType row = *groupStartIt; ValueType xBest, yBest; - multiplyRow(row, A, multiplier, b[row], xBest, yBest); + multiplyRow(row, b[row], xBest, yBest); ++row; // Only do more work if there are still rows in this row group if (row != groupEnd) { @@ -697,7 +727,7 @@ namespace storm { ValueType bestValue = xBest + yBest * getPrimaryBound(); for (;row < groupEnd; ++row) { // Get the multiplication results - multiplyRow(row, A, multiplier, b[row], xi, yi); + multiplyRow(row, b[row], xi, yi); ValueType currentValue = xi + yi * getPrimaryBound(); // Check if the current row is better then the previously found one if (better(currentValue, bestValue)) { @@ -718,17 +748,17 @@ namespace storm { } template - void performIterationStepUpdateDecisionValue(storm::storage::SparseMatrix const& A, storm::solver::Multiplier const& multiplier, std::vector const& b) { + void performIterationStepUpdateDecisionValue(std::vector const& b) { auto xIt = x.rbegin(); auto yIt = y.rbegin(); - auto groupStartIt = A.getRowGroupIndices().rbegin(); + auto groupStartIt = rowGroupIndices.rbegin(); uint64_t groupEnd = *groupStartIt; ++groupStartIt; - for (auto groupStartIte = A.getRowGroupIndices().rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { + for (auto groupStartIte = rowGroupIndices.rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { // Perform the iteration for the first row in the group uint64_t row = *groupStartIt; ValueType xBest, yBest; - multiplyRow(row, A, multiplier, b[row], xBest, yBest); + multiplyRow(row, b[row], xBest, yBest); ++row; // Only do more work if there are still rows in this row group if (row != groupEnd) { @@ -738,7 +768,7 @@ namespace storm { ValueType bestValue = xBest + yBest * getPrimaryBound(); for (;row < groupEnd; ++row) { // Get the multiplication results - multiplyRow(row, A, multiplier, b[row], xi, yi); + multiplyRow(row, b[row], xi, yi); ValueType currentValue = xi + yi * getPrimaryBound(); // Check if the current row is better then the previously found one if (better(currentValue, bestValue)) { @@ -765,7 +795,7 @@ namespace storm { } } else { for (;row < groupEnd; ++row) { - multiplyRow(row, A, multiplier, b[row], xi, yi); + multiplyRow(row, b[row], xi, yi); // Update the best choice if (yi > yBest || (yi == yBest && better(xi, xBest))) { xTmp[xyTmpIndex] = std::move(xBest); @@ -963,6 +993,12 @@ namespace storm { bool relative; ValueType precision; + + std::vector matrixValues; + std::vector matrixColumns; + std::vector rowIndications; + std::vector::index_type> const& rowGroupIndices; + IndexType numRows; }; template @@ -973,12 +1009,9 @@ namespace storm { if (!this->auxiliaryRowGroupVector) { this->auxiliaryRowGroupVector = std::make_unique>(); } - - if (!this->multiplierA) { - this->multiplierA = storm::solver::MultiplierFactory().create(env, *this->A); - } - - SoundValueIterationHelper helper(x, *this->auxiliaryRowGroupVector, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().minMax().getPrecision()), this->A->getSizeOfLargestRowGroup()); + + // TODO: implement caching for the helper + SoundValueIterationHelper helper(*this->A, x, *this->auxiliaryRowGroupVector, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().minMax().getPrecision())); // Prepare initial bounds for the solution (if given) if (this->hasLowerBound()) { @@ -999,13 +1032,13 @@ namespace storm { while (status == SolverStatus::InProgress && iterations < env.solver().minMax().getMaximalNumberOfIterations()) { if (minimize(dir)) { - helper.template performIterationStep(*this->A, *this->multiplierA, b); + helper.template performIterationStep(b); if (helper.template checkConvergenceUpdateBounds(relevantValuesPtr)) { status = SolverStatus::Converged; } } else { assert(maximize(dir)); - helper.template performIterationStep(*this->A, *this->multiplierA, b); + helper.template performIterationStep(b); if (helper.template checkConvergenceUpdateBounds(relevantValuesPtr)) { status = SolverStatus::Converged; } diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 48a2819c7..96daa55b4 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -1,5 +1,7 @@ #include "storm/solver/NativeLinearEquationSolver.h" +#include + #include "storm/environment/solver/NativeSolverEnvironment.h" #include "storm/utility/ConstantsComparator.h" @@ -12,7 +14,7 @@ #include "storm/exceptions/InvalidEnvironmentException.h" #include "storm/exceptions/UnmetRequirementException.h" #include "storm/exceptions/PrecisionExceededException.h" - +#include "storm/exceptions/NotSupportedException.h" namespace storm { namespace solver { @@ -573,11 +575,31 @@ namespace storm { template class SoundPowerHelper { public: - SoundPowerHelper(std::vector& x, std::vector& y, bool relative, ValueType const& precision) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), minIndex(0), maxIndex(0), relative(relative), precision(precision) { + + typedef uint32_t IndexType; + + SoundPowerHelper(storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector& y, bool relative, ValueType const& precision) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), minIndex(0), maxIndex(0), relative(relative), precision(precision) { + STORM_LOG_THROW(matrix.getEntryCount() < std::numeric_limits::max(), storm::exceptions::NotSupportedException, "The number of matrix entries is too large for the selected index type."); x.assign(x.size(), storm::utility::zero()); y.assign(x.size(), storm::utility::one()); convergencePhase1 = true; firstIndexViolatingConvergence = 0; + + numRows = matrix.getRowCount(); + matrixValues.clear(); + matrixColumns.clear(); + rowIndications.clear(); + matrixValues.reserve(matrix.getNonzeroEntryCount()); + matrixColumns.reserve(matrix.getColumnCount()); + rowIndications.reserve(numRows + 1); + rowIndications.push_back(0); + for (IndexType r = 0; r < numRows; ++r) { + for (auto const& entry : matrix.getRow(r)) { + matrixValues.push_back(entry.getValue()); + matrixColumns.push_back(entry.getColumn()); + } + rowIndications.push_back(matrixValues.size()); + } } inline void setLowerBound(ValueType const& value) { @@ -590,21 +612,29 @@ namespace storm { upperBound = value; } - void multiplyRow(uint64_t const& row, storm::storage::SparseMatrix const& A, storm::solver::Multiplier const& multiplier, ValueType const& bi, ValueType& xi, ValueType& yi) { + void multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi) { + assert(rowIndex < numRows); ValueType xRes = bi; ValueType yRes = storm::utility::zero(); - multiplier.multiplyRow2(row, x, xRes, y, yRes); + + auto entryIt = matrixValues.begin() + rowIndications[rowIndex]; + auto entryItE = matrixValues.begin() + rowIndications[rowIndex + 1]; + auto colIt = matrixColumns.begin() + rowIndications[rowIndex]; + for (; entryIt != entryItE; ++entryIt, ++colIt) { + xRes += *entryIt * x[*colIt]; + yRes += *entryIt * y[*colIt]; + } xi = std::move(xRes); yi = std::move(yRes); } - void performIterationStep(storm::storage::SparseMatrix const& A, storm::solver::Multiplier const& multiplier, std::vector const& b) { + void performIterationStep(std::vector const& b) { auto xIt = x.rbegin(); auto yIt = y.rbegin(); - uint64_t row = A.getRowCount(); + IndexType row = numRows; while (row > 0) { --row; - multiplyRow(row, A, multiplier, b[row], *xIt, *yIt); + multiplyRow(row, b[row], *xIt, *yIt); ++xIt; ++yIt; } @@ -744,6 +774,11 @@ namespace storm { bool convergencePhase1; uint64_t firstIndexViolatingConvergence; + std::vector matrixValues; + std::vector matrixColumns; + std::vector rowIndications; + IndexType numRows; + bool relative; ValueType precision; }; @@ -757,11 +792,8 @@ namespace storm { this->cachedRowVector = std::make_unique>(); } - if (!this->multiplier) { - this->multiplier = storm::solver::MultiplierFactory().create(env, *this->A); - } - - SoundPowerHelper helper(x, *this->cachedRowVector, env.solver().native().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().native().getPrecision())); + // TODO: implement caching for the helper + SoundPowerHelper helper(*this->A, x, *this->cachedRowVector, env.solver().native().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().native().getPrecision())); // Prepare initial bounds for the solution (if given) if (this->hasLowerBound()) { @@ -782,7 +814,7 @@ namespace storm { uint64_t iterations = 0; while (!converged && iterations < env.solver().native().getMaximalNumberOfIterations()) { - helper.performIterationStep(*this->A, *this->multiplier, b); + helper.performIterationStep(b); if (helper.checkConvergenceUpdateBounds(relevantValuesPtr)) { converged = true; }