Multiplier: Added a flag to specify whether gaussSeidel style multiplications should be performed forward or backwards.

src/storm/solver/AcyclicLinearEquationSolver.h

@@ -0,0 +1,47 @@
+#pragma once
+
+#include "storm/solver/StandardMinMaxLinearEquationSolver.h"
+
+namespace storm {
+
+    class Environment;
+
+    namespace solver {
+
+        /*!
+         * This solver can be used on equation systems that are known to be acyclic.
+         * It is optimized for solving many instances of the equation system with the same underlying matrix.
+         */
+        template<typename ValueType>
+        class AcyclicMinMaxLinearEquationSolver : public StandardMinMaxLinearEquationSolver<ValueType> {
+        public:
+            AcyclicMinMaxLinearEquationSolver();
+            AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
+            AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A);
+
+            virtual ~AcyclicMinMaxLinearEquationSolver() {
+            }
+
+            virtual void clearCache() const override;
+
+            virtual MinMaxLinearEquationSolverRequirements getRequirements(Environment const& env, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none, bool const& hasInitialScheduler = false) const override ;
+
+        protected:
+
+            virtual bool internalSolveEquations(storm::Environment const& env, OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
+
+        private:
+            // cached multiplier either with original matrix or ordered matrix
+            mutable std::unique_ptr<storm::solver::Multiplier<ValueType>> multiplier;
+            // cached matrix for the multiplier (only if different from original matrix)
+            mutable boost::optional<storm::storage::SparseMatrix<ValueType>> orderedMatrix;
+            // cached row group ordering (only if not identity)
+            mutable boost::optional<std::vector<uint64_t>> rowGroupOrdering; // A.rowGroupCount() entries
+            // can be used if the entries in 'b' need to be reordered
+            mutable boost::optional<std::vector<ValueType>> auxiliaryRowVector; // A.rowCount() entries
+            // contains factors applied to scale the entries of the 'b' vector
+            mutable std::vector<std::pair<uint64_t, ValueType>> bFactors;
+
+        };
+    }
+}

src/storm/solver/GmmxxMultiplier.cpp

@@ -65,14 +65,22 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void GmmxxMultiplier<ValueType>::multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
+        void GmmxxMultiplier<ValueType>::multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, bool backwards) const {
             initialize();
             STORM_LOG_ASSERT(gmmMatrix.nr == gmmMatrix.nc, "Expecting square matrix.");
+            if (backwards) {
                 if (b) {
                     gmm::mult_add_by_row_bwd(gmmMatrix, x, *b, x, gmm::abstract_dense());
                 } else {
                     gmm::mult_by_row_bwd(gmmMatrix, x, x, gmm::abstract_dense());
                 }
+            } else {
+                if (b) {
+                    gmm::mult_add_by_row(gmmMatrix, x, *b, x, gmm::abstract_dense());
+                } else {
+                    gmm::mult_by_row(gmmMatrix, x, x, gmm::abstract_dense());
+                }
+            }
         }
         
         template<typename ValueType>
@@ -90,7 +98,7 @@ namespace storm {
             if (parallelize(env)) {
                 multAddReduceParallel(dir, rowGroupIndices, x, b, *target, choices);
             } else {
-                multAddReduceHelper(dir, rowGroupIndices, x, b, *target, choices);
+                multAddReduceHelper(dir, rowGroupIndices, x, b, *target, choices, false);
             }
             if (&x == &result) {
                 std::swap(result, *this->cachedVector);
@@ -98,9 +106,9 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void GmmxxMultiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices) const {
+        void GmmxxMultiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices, bool backwards) const {
             initialize();
-            multAddReduceHelper(dir, rowGroupIndices, x, b, x, choices);
+            multAddReduceHelper(dir, rowGroupIndices, x, b, x, choices, backwards);
         }
         
         template<typename ValueType>
@@ -119,16 +127,24 @@ namespace storm {
         }
 
         template<typename ValueType>
-        void GmmxxMultiplier<ValueType>::multAddReduceHelper(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices) const {
+        void GmmxxMultiplier<ValueType>::multAddReduceHelper(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices, bool backwards) const {
             if (dir == storm::OptimizationDirection::Minimize) {
-                multAddReduceHelper<storm::utility::ElementLess<ValueType>>(rowGroupIndices, x, b, result, choices);
+                if (backwards) {
+                    multAddReduceHelper<storm::utility::ElementLess<ValueType>, true>(rowGroupIndices, x, b, result, choices);
                 } else {
-                multAddReduceHelper<storm::utility::ElementGreater<ValueType>>(rowGroupIndices, x, b, result, choices);
+                    multAddReduceHelper<storm::utility::ElementLess<ValueType>, false>(rowGroupIndices, x, b, result, choices);
+                }
+            } else {
+                if (backwards) {
+                    multAddReduceHelper<storm::utility::ElementGreater<ValueType>, true>(rowGroupIndices, x, b, result, choices);
+                } else {
+                    multAddReduceHelper<storm::utility::ElementGreater<ValueType>, false>(rowGroupIndices, x, b, result, choices);
+                }
             }
         }
         
         template<typename ValueType>
-        template<typename Compare>
+        template<typename Compare, bool backwards>
         void GmmxxMultiplier<ValueType>::multAddReduceHelper(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices) const {
             Compare compare;
             typedef std::vector<ValueType> VectorType;
@@ -136,29 +152,31 @@ namespace storm {
             
             typename gmm::linalg_traits<VectorType>::const_iterator add_it, add_ite;
             if (b) {
-                add_it = gmm::vect_end(*b) - 1;
-                add_ite = gmm::vect_begin(*b) - 1;
+                add_it = backwards ? gmm::vect_end(*b) - 1 : gmm::vect_begin(*b);
+                add_ite = backwards ?gmm::vect_begin(*b) - 1 : gmm::vect_end(*b);
             }
-            typename gmm::linalg_traits<VectorType>::iterator target_it = gmm::vect_end(result) - 1;
-            typename gmm::linalg_traits<MatrixType>::const_row_iterator itr = mat_row_const_end(gmmMatrix) - 1;
+            typename gmm::linalg_traits<VectorType>::iterator target_it = backwards ? gmm::vect_end(result) - 1 : gmm::vect_begin(result);
+            typename gmm::linalg_traits<MatrixType>::const_row_iterator itr = backwards ? mat_row_const_end(gmmMatrix) - 1 : mat_row_const_begin(gmmMatrix);
             typename std::vector<uint64_t>::iterator choice_it;
             if (choices) {
-                choice_it = choices->end() - 1;
+                choice_it = backwards ? choices->end() - 1 : choices->begin();
             }
 
             // Variables for correctly tracking choices (only update if new choice is strictly better).
             ValueType oldSelectedChoiceValue;
             uint64_t selectedChoice;
 
-            uint64_t currentRow = gmmMatrix.nrows() - 1;
-            for (auto row_group_it = rowGroupIndices.end() - 2, row_group_ite = rowGroupIndices.begin() - 1; row_group_it != row_group_ite; --row_group_it, --choice_it, --target_it) {
+            uint64_t currentRow = backwards ? gmmMatrix.nrows() - 1 : 0;
+            auto row_group_it = backwards ? rowGroupIndices.end() - 2 : rowGroupIndices.begin();
+            auto row_group_ite = backwards ? rowGroupIndices.begin() - 1 : rowGroupIndices.end() - 1;
+            while (row_group_it != row_group_ite) {
                 ValueType currentValue = storm::utility::zero<ValueType>();
                 
                 // Only multiply and reduce if the row group is not empty.
                 if (*row_group_it != *(row_group_it + 1)) {
+                    // Process the (backwards ? last : first) row of the current row group
                     if (b) {
                         currentValue = *add_it;
-                        --add_it;
                     }
                     
                     currentValue += vect_sp(gmm::linalg_traits<MatrixType>::row(itr), x);
@@ -170,10 +188,20 @@ namespace storm {
                         }
                     }
                     
+                    // move row-based iterators to the next row
+                    if (backwards) {
                         --itr;
                         --currentRow;
+                        --add_it;
+                    } else {
+                        ++itr;
+                        ++currentRow;
+                        ++add_it;
+                    }
                     
-                    for (uint64_t row = *row_group_it + 1, rowEnd = *(row_group_it + 1); row < rowEnd; ++row, --currentRow, --itr, --add_it) {
+                    // Process the (rowGroupSize-1) remaining rows within the current row Group
+                    uint64_t rowGroupSize = *(row_group_it + 1) - *row_group_it;
+                    for (uint64_t i = 1; i < rowGroupSize; ++i) {
                         ValueType newValue = b ? *add_it : storm::utility::zero<ValueType>();
                         newValue += vect_sp(gmm::linalg_traits<MatrixType>::row(itr), x);
                         
@@ -187,6 +215,16 @@ namespace storm {
                                 selectedChoice = currentRow - *row_group_it;
                             }
                         }
+                        // move row-based iterators to the next row
+                        if (backwards) {
+                             --itr;
+                             --currentRow;
+                             --add_it;
+                        } else {
+                             ++itr;
+                             ++currentRow;
+                             ++add_it;
+                        }
                     }
                     
                     // Finally write value to target vector.
@@ -195,11 +233,22 @@ namespace storm {
                         *choice_it = selectedChoice;
                     }
                 }
+                
+                // move rowGroup-based iterators to the next row group
+                if (backwards) {
+                    --row_group_it;
+                    --choice_it;
+                    --target_it;
+                } else {
+                    ++row_group_it;
+                    ++choice_it;
+                    ++target_it;
+                }
             }
         }
         
         template<>
-        template<typename Compare>
+        template<typename Compare, bool backwards>
         void GmmxxMultiplier<storm::RationalFunction>::multAddReduceHelper(std::vector<uint64_t> const& rowGroupIndices, std::vector<storm::RationalFunction> const& x, std::vector<storm::RationalFunction> const* b, std::vector<storm::RationalFunction>& result, std::vector<uint64_t>* choices) const {
             STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Operation not supported for this data type.");
         }

src/storm/solver/GmmxxMultiplier.h

@@ -22,9 +22,9 @@ namespace storm {
             virtual ~GmmxxMultiplier() = default;
             
             virtual void multiply(Environment const& env, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
-            virtual void multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b) const override;
+            virtual void multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, bool backwards = true) const override;
             virtual void multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices = nullptr) const override;
-            virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr) const override;
+            virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr, bool backwards = true) const override;
             virtual void multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x, ValueType& value) const override;
             virtual void clearCache() const override;
             
@@ -36,9 +36,9 @@ namespace storm {
             void multAdd(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const;
             void multAddParallel(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const;
             void multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices = nullptr) const;
-            void multAddReduceHelper(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices = nullptr) const;
+            void multAddReduceHelper(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices = nullptr, bool backwards = true) const;
             
-            template<typename Compare>
+            template<typename Compare, bool backwards = true>
             void multAddReduceHelper(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices = nullptr) const;
 
             mutable gmm::csr_matrix<ValueType> gmmMatrix;

src/storm/solver/Multiplier.cpp

@@ -33,8 +33,8 @@ namespace storm {
         }
 
         template<typename ValueType>
-        void Multiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices) const {
-            multiplyAndReduceGaussSeidel(env, dir, this->matrix.getRowGroupIndices(), x, b, choices);
+        void Multiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices, bool backwards) const {
+            multiplyAndReduceGaussSeidel(env, dir, this->matrix.getRowGroupIndices(), x, b, choices, backwards);
         }
     
         template<typename ValueType>

src/storm/solver/Multiplier.h

@@ -49,8 +49,9 @@ namespace storm {
              * to the number of columns of A.
              * @param b If non-null, this vector is added after the multiplication. If given, its length must be equal
              * to the number of rows of A.
+             * @param backwards if true, the iterations will be performed beginning from the last row and ending at the first row.
              */
-            virtual void multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b) const = 0;
+            virtual void multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, bool backwards = true) const = 0;
             
             /*!
              * Performs a matrix-vector multiplication x' = A*x + b and then minimizes/maximizes over the row groups
@@ -82,9 +83,10 @@ namespace storm {
              * @param result The target vector into which to write the multiplication result. Its length must be equal
              * to the number of rows of A. Can be the same as the x vector.
              * @param choices If given, the choices made in the reduction process are written to this vector.
+             * @param backwards if true, the iterations will be performed beginning from the last rowgroup and ending at the first rowgroup.
              */
-            void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr) const;
-            virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr) const = 0;
+            void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr, bool backwards = true) const;
+            virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr, bool backwards = true) const = 0;
             
             /*!
              * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After

src/storm/solver/NativeMultiplier.cpp

@@ -53,8 +53,12 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
+        void NativeMultiplier<ValueType>::multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, bool backwards) const {
+            if (backwards) {
                 this->matrix.multiplyWithVectorBackward(x, x, b);
+            } else {
+                this->matrix.multiplyWithVectorForward(x, x, b);
+            }
         }
         
         template<typename ValueType>
@@ -79,8 +83,12 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices) const {
+        void NativeMultiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices, bool backwards) const {
+            if (backwards) {
                 this->matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices);
+            } else {
+                this->matrix.multiplyAndReduceForward(dir, rowGroupIndices, x, b, x, choices);
+            }
         }
         
         template<typename ValueType>

src/storm/solver/NativeMultiplier.h

@@ -19,9 +19,9 @@ namespace storm {
             virtual ~NativeMultiplier() = default;
             
             virtual void multiply(Environment const& env, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
-            virtual void multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b) const override;
+            virtual void multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, bool backwards = true) const override;
             virtual void multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices = nullptr) const override;
-            virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr) const override;
+            virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr, bool backwards = true) const override;
             virtual void multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x, ValueType& value) const override;
             virtual void multiplyRow2(uint64_t const& rowIndex, std::vector<ValueType> const& x1, ValueType& val1, std::vector<ValueType> const& x2, ValueType& val2) const override;