trying an alternative implementation of qvi

src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

@@ -679,7 +679,7 @@ namespace storm {
             template<OptimizationDirection dir>
             void performIterationStep(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType> const& b) {
                 if (!decisionValueBlocks) {
-                    performIterationStepUpdateDecisionValue<dir>(A, b);
+                    performIterationStepUpdateDecisionValue2<dir>(A, b);
                 } else {
                     assert(decisionValue == getPrimaryBound<dir>());
                     auto xIt = x.rbegin();
@@ -800,6 +800,86 @@ namespace storm {
                 }
             }
             
+            template<OptimizationDirection dir>
+            void performIterationStepUpdateDecisionValue2(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType> const& b) {
+                auto const& groupIndices = A.getRowGroupIndices();
+                uint64_t group = groupIndices.size();
+                while (group != 0) {
+                    uint64_t const& groupEnd = groupIndices[group];
+                    --group;
+                    uint64_t row = groupIndices[group];
+                    ValueType xBest, yBest;
+                    multiplyRow(row, A, b[row], xBest, yBest);
+                    ++row;
+                    
+                    // Only do more work if there are still rows in this row group
+                    if (row != groupEnd) {
+                        ValueType xi, yi;
+                        uint64_t xyTmpIndex = 0;
+                        if (hasPrimaryBound<dir>()) {
+                            ValueType bestValue = xBest + yBest * getPrimaryBound<dir>();
+                            for (;row < groupEnd; ++row) {
+                                // Get the multiplication results
+                                multiplyRow(row, A, b[row], xi, yi);
+                                ValueType currentValue = xi + yi * getPrimaryBound<dir>();
+                                // Check if the current row is better then the previously found one
+                                if (better<dir>(currentValue, bestValue)) {
+                                    if (yBest < yi) {
+                                        // We need to store the 'old' best value as it might be relevant for the decision value
+                                        xTmp[xyTmpIndex] = std::move(xBest);
+                                        yTmp[xyTmpIndex] = std::move(yBest);
+                                        ++xyTmpIndex;
+                                    }
+                                    xBest = std::move(xi);
+                                    yBest = std::move(yi);
+                                    bestValue = std::move(currentValue);
+                                } else if (yBest > yi) {
+                                    // If the value for this row is not strictly better, it might still be equal and have a better y value
+                                    if (currentValue == bestValue) {
+                                        xBest = std::move(xi);
+                                        yBest = std::move(yi);
+                                    } else {
+                                        xTmp[xyTmpIndex] = std::move(xi);
+                                        yTmp[xyTmpIndex] = std::move(yi);
+                                        ++xyTmpIndex;
+                                    }
+                                }
+                            }
+                        } else {
+                            for (;row < groupEnd; ++row) {
+                                multiplyRow(row, A, b[row], xi, yi);
+                                // Update the best choice
+                                if (yi > yBest || (yi == yBest && better<dir>(xi, xBest))) {
+                                        xTmp[xyTmpIndex] = std::move(xBest);
+                                        yTmp[xyTmpIndex] = std::move(yBest);
+                                        ++xyTmpIndex;
+                                    xBest = std::move(xi);
+                                    yBest = std::move(yi);
+                                } else {
+                                    xTmp[xyTmpIndex] = std::move(xi);
+                                    yTmp[xyTmpIndex] = std::move(yi);
+                                    ++xyTmpIndex;
+                                }
+                            }
+                        }
+                        
+                        // Update the decision value
+                        for (uint64_t i = 0; i < xyTmpIndex; ++i) {
+                            ValueType deltaY = yBest - yTmp[i];
+                            if (deltaY > storm::utility::zero<ValueType>()) {
+                                ValueType newDecisionValue = (xTmp[i] - xBest) / deltaY;
+                                if (!hasDecisionValue || better<dir>(newDecisionValue, decisionValue)) {
+                                    decisionValue = std::move(newDecisionValue);
+                                    hasDecisionValue = true;
+                                }
+                            }
+                        }
+                    }
+                    x[group] = std::move(xBest);
+                    y[group] = std::move(yBest);
+                }
+            }
+            
             bool checkRestartCriterion() {
                 return false;
                 // iterations <= restartMaxIterations && (minimize(dir) ? restartThreshold * improvedPrimaryBound > primaryBound : restartThreshold * primaryBound > improvedPrimaryBound