allow for summand in matrix-vector multiplication

src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

@@ -256,7 +256,6 @@ namespace storm {
                 storm::utility::vector::selectVectorValues<ValueType>(*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()); }

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -207,22 +207,15 @@ namespace storm {
         template<typename ValueType>
         void NativeLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& 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) {
                     this->cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount());
                 }
                 
-                A->multiplyWithVector(x, *this->cachedRowVector);
-                if (b != nullptr) {
-                    storm::utility::vector::addVectors(*this->cachedRowVector, *b, result);
-                } else {
+                A->multiplyWithVector(x, *this->cachedRowVector, b);
                 result.swap(*this->cachedRowVector);
-                }
                 
                 if (!this->isCachingEnabled()) {
                     clearCache();

src/storm/storage/SparseMatrix.cpp

@@ -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;
         }
         
@@ -1293,20 +1291,20 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void SparseMatrix<ValueType>::multiplyWithVector(std::vector<ValueType> const& vector, std::vector<ValueType>& result) const {
+        void SparseMatrix<ValueType>::multiplyWithVector(std::vector<ValueType> const& vector, std::vector<ValueType>& result, std::vector<value_type> 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<typename ValueType>
-        void SparseMatrix<ValueType>::multiplyWithVectorSequential(std::vector<ValueType> const& vector, std::vector<ValueType>& result) const {
+        void SparseMatrix<ValueType>::multiplyWithVectorSequential(std::vector<ValueType> const& vector, std::vector<ValueType>& result, std::vector<value_type> 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<ValueType> tmpVector(this->getRowCount());
@@ -1318,9 +1316,18 @@ namespace storm {
                 std::vector<index_type>::const_iterator rowIterator = rowIndications.begin();
                 typename std::vector<ValueType>::iterator resultIterator = result.begin();
                 typename std::vector<ValueType>::iterator resultIteratorEnd = result.end();
+                typename std::vector<ValueType>::const_iterator summandIterator;
+                if (summand) {
+                    summandIterator = summand->begin();
+                }
                 
                 for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) {
+                    if (summand) {
+                        *resultIterator = *summandIterator;
+                        ++summandIterator;
+                    } else {
                         *resultIterator = storm::utility::zero<ValueType>();
+                    }
                     
                     for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
                         *resultIterator += it->getValue() * vector[it->getColumn()];
@@ -1331,7 +1338,7 @@ namespace storm {
         
 #ifdef STORM_HAVE_INTELTBB
         template<typename ValueType>
-        void SparseMatrix<ValueType>::multiplyWithVectorParallel(std::vector<ValueType> const& vector, std::vector<ValueType>& result) const {
+        void SparseMatrix<ValueType>::multiplyWithVectorParallel(std::vector<ValueType> const& vector, std::vector<ValueType>& result, std::vector<value_type> 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<ValueType> tmpVector(this->getRowCount());
@@ -1348,9 +1355,18 @@ namespace storm {
                                       std::vector<index_type>::const_iterator rowIteratorEnd = this->rowIndications.begin() + endRow;
                                       typename std::vector<ValueType>::iterator resultIterator = result.begin() + startRow;
                                       typename std::vector<ValueType>::iterator resultIteratorEnd = result.begin() + endRow;
+                                      typename std::vector<ValueType>::const_iterator summandIterator;
+                                      if (summand) {
+                                          summandIterator = summand->begin() + startRow;
+                                      }
                                       
                                       for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) {
+                                          if (summand) {
+                                              *resultIterator = *summandIterator;
+                                              ++summandIterator;
+                                          } else {
                                               *resultIterator = storm::utility::zero<ValueType>();
+                                          }
                                           
                                           for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
                                               *resultIterator += it->getValue() * vector[it->getColumn()];

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<value_type> const& vector, std::vector<value_type>& result) const;
+            void multiplyWithVector(std::vector<value_type> const& vector, std::vector<value_type>& result, std::vector<value_type> 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<value_type> const& vector, std::vector<value_type>& result) const;
+            void multiplyWithVectorSequential(std::vector<value_type> const& vector, std::vector<value_type>& result, std::vector<value_type> 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<value_type> const& vector, std::vector<value_type>& result) const;
+            void multiplyWithVectorParallel(std::vector<value_type> const& vector, std::vector<value_type>& result, std::vector<value_type> 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<std::vector<index_type>> rowGroupIndices;
-            
         };
         
 #ifdef STORM_HAVE_CARL

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<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
     
     ASSERT_NO_THROW(solver->solveEquations(storm::OptimizationDirection::Maximize, x, b));