started to fill value iteration implementation in new general min-max solver

Former-commit-id: e54cb8a0f9

src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp

@@ -39,7 +39,6 @@ namespace storm {
                 maybeStates &= ~psiStates;
                 STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
                 
-                
                 if (!maybeStates.empty()) {
                     // We can eliminate the rows and columns from the original transition probability matrix that have probability 0.
                     storm::storage::SparseMatrix<ValueType> submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, false);

src/settings/SettingsManager.cpp

@@ -23,11 +23,12 @@
 #include "src/settings/modules/EigenEquationSolverSettings.h"
 #include "src/settings/modules/GmmxxEquationSolverSettings.h"
 #include "src/settings/modules/NativeEquationSolverSettings.h"
+#include "src/settings/modules/EliminationSettings.h"
+#include "src/settings/modules/MinMaxEquationSolverSettings.h"
 #include "src/settings/modules/BisimulationSettings.h"
 #include "src/settings/modules/GlpkSettings.h"
 #include "src/settings/modules/GurobiSettings.h"
 #include "src/settings/modules/ParametricSettings.h"
-#include "src/settings/modules/EliminationSettings.h"
 #include "src/settings/modules/TopologicalValueIterationEquationSolverSettings.h"
 #include "src/settings/modules/ExplorationSettings.h"
 #include "src/utility/macros.h"
@@ -511,12 +512,13 @@ namespace storm {
             storm::settings::addModule<storm::settings::modules::GmmxxEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::EigenEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::NativeEquationSolverSettings>();
+            storm::settings::addModule<storm::settings::modules::EliminationSettings>();
+            storm::settings::addModule<storm::settings::modules::MinMaxEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::BisimulationSettings>();
             storm::settings::addModule<storm::settings::modules::GlpkSettings>();
             storm::settings::addModule<storm::settings::modules::GurobiSettings>();
             storm::settings::addModule<storm::settings::modules::TopologicalValueIterationEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::ParametricSettings>();
-            storm::settings::addModule<storm::settings::modules::EliminationSettings>();
             storm::settings::addModule<storm::settings::modules::ExplorationSettings>();
         }
 

src/settings/modules/CoreSettings.cpp

@@ -142,20 +142,6 @@ namespace storm {
                 this->engine = newEngine;
             }
             
-            storm::solver::MinMaxTechnique CoreSettings::getMinMaxEquationSolvingTechnique() const {
-				std::string minMaxEquationSolvingTechnique = this->getOption(minMaxEquationSolvingTechniqueOptionName).getArgumentByName("name").getValueAsString();
-				if (minMaxEquationSolvingTechnique == "value-iteration" || minMaxEquationSolvingTechnique == "vi") {
-					return storm::solver::MinMaxTechnique::ValueIteration;
-				} else if (minMaxEquationSolvingTechnique == "policy-iteration" || minMaxEquationSolvingTechnique == "pi") {
-					return storm::solver::MinMaxTechnique::PolicyIteration;
-				}
-				STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown min/max equation solving technique '" << minMaxEquationSolvingTechnique << "'.");
-			}
-
-			bool CoreSettings::isMinMaxEquationSolvingTechniqueSet() const {
-				return this->getOption(minMaxEquationSolvingTechniqueOptionName).getHasOptionBeenSet();
-			}
-            
             void CoreSettings::finalize() {
                 // Finalize engine.
                 std::string engineStr = this->getOption(engineOptionName).getArgumentByName("name").getValueAsString();

src/settings/modules/CoreSettings.h

@@ -10,7 +10,7 @@ namespace storm {
     namespace solver {
         enum class EquationSolverType;
         enum class LpSolverType;
-        enum class MinMaxTechnique;
+        enum class MinMaxMethod;
         enum class SmtSolverType;
     }
     
@@ -22,7 +22,7 @@ namespace storm {
         namespace modules {
 
             /*!
-             * This class represents the markov chain settings.
+             * This class represents the core settings.
              */
             class CoreSettings : public ModuleSettings {
             public:
@@ -32,7 +32,7 @@ namespace storm {
                 };
 
                 /*!
-                 * Creates a new set of markov chain settings.
+                 * Creates a new set of core settings.
                  */
                 CoreSettings();
 
@@ -128,21 +128,6 @@ namespace storm {
                  */
                 void setEngine(Engine);
 
-                /*!
-                 * Retrieves whether a min/max equation solving technique has been set.
-                 *
-                 * @return True iff an equation solving technique has been set.
-                 */
-                bool isMinMaxEquationSolvingTechniqueSet() const;
-
-                /*!
-                 * Retrieves the selected min/max equation solving technique.
-                 *
-                 * @return The selected min/max equation solving technique.
-                 */
-                storm::solver::MinMaxTechnique getMinMaxEquationSolvingTechnique() const;
-
-
                 bool check() const override;
                 void finalize() override;
 

src/settings/modules/MinMaxEquationSolverSettings.cpp

@@ -0,0 +1,72 @@
+#include "src/settings/modules/MinMaxEquationSolverSettings.h"
+
+#include "src/settings/Option.h"
+#include "src/settings/ArgumentBuilder.h"
+#include "src/settings/OptionBuilder.h"
+
+#include "src/utility/macros.h"
+#include "src/exceptions/IllegalArgumentValueException.h"
+
+namespace storm {
+    namespace settings {
+        namespace modules {
+            
+            const std::string MinMaxEquationSolverSettings::moduleName = "minmax";
+            const std::string MinMaxEquationSolverSettings::solvingMethodOptionName = "method";
+            const std::string MinMaxEquationSolverSettings::maximalIterationsOptionName = "maxiter";
+            const std::string MinMaxEquationSolverSettings::maximalIterationsOptionShortName = "i";
+            const std::string MinMaxEquationSolverSettings::precisionOptionName = "precision";
+
+            MinMaxEquationSolverSettings::MinMaxEquationSolverSettings() : ModuleSettings(moduleName) {
+                std::vector<std::string> minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration"};
+                this->addOption(storm::settings::OptionBuilder(moduleName, solvingMethodOptionName, false, "Sets which min/max linear equation solving technique is preferred.")
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a min/max linear equation solving technique. Available are: value-iteration (vi) and policy-iteration (pi).").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(minMaxSolvingTechniques)).setDefaultValueString("vi").build()).build());
+                
+                this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build());
+                
+                this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+
+                this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build());
+            }
+            
+            storm::solver::MinMaxMethod MinMaxEquationSolverSettings::getMinMaxEquationSolvingMethod() const {
+                std::string minMaxEquationSolvingTechnique = this->getOption(solvingMethodOptionName).getArgumentByName("name").getValueAsString();
+                if (minMaxEquationSolvingTechnique == "value-iteration" || minMaxEquationSolvingTechnique == "vi") {
+                    return storm::solver::MinMaxMethod::ValueIteration;
+                } else if (minMaxEquationSolvingTechnique == "policy-iteration" || minMaxEquationSolvingTechnique == "pi") {
+                    return storm::solver::MinMaxMethod::PolicyIteration;
+                }
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown min/max equation solving technique '" << minMaxEquationSolvingTechnique << "'.");
+            }
+            
+            bool MinMaxEquationSolverSettings::isMinMaxEquationSolvingMethodSet() const {
+                return this->getOption(solvingMethodOptionName).getHasOptionBeenSet();
+            }
+            
+            bool MinMaxEquationSolverSettings::isMaximalIterationCountSet() const {
+                return this->getOption(maximalIterationsOptionName).getHasOptionBeenSet();
+            }
+            
+            uint_fast64_t MinMaxEquationSolverSettings::getMaximalIterationCount() const {
+                return this->getOption(maximalIterationsOptionName).getArgumentByName("count").getValueAsUnsignedInteger();
+            }
+            
+            bool MinMaxEquationSolverSettings::isPrecisionSet() const {
+                return this->getOption(precisionOptionName).getHasOptionBeenSet();
+            }
+            
+            double MinMaxEquationSolverSettings::getPrecision() const {
+                return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble();
+            }
+            
+            bool MinMaxEquationSolverSettings::isConvergenceCriterionSet() const {
+                return this->getOption(absoluteOptionName).getHasOptionBeenSet();
+            }
+            
+            MinMaxEquationSolverSettings::ConvergenceCriterion MinMaxEquationSolverSettings::getConvergenceCriterion() const {
+                return this->getOption(absoluteOptionName).getHasOptionBeenSet() ? MinMaxEquationSolverSettings::ConvergenceCriterion::Absolute : MinMaxEquationSolverSettings::ConvergenceCriterion::Relative;
+            }
+            
+        }
+    }
+}

src/settings/modules/MinMaxEquationSolverSettings.h

@@ -0,0 +1,91 @@
+#pragma once
+
+#include "storm-config.h"
+#include "src/settings/modules/ModuleSettings.h"
+
+#include "src/solver/SolverSelectionOptions.h"
+
+namespace storm {
+    namespace settings {
+        namespace modules {
+            
+            /*!
+             * This class represents the min/max solver settings.
+             */
+            class MinMaxEquationSolverSettings : public ModuleSettings {
+            public:
+                // An enumeration of all available convergence criteria.
+                enum class ConvergenceCriterion { Absolute, Relative };
+                
+                MinMaxEquationSolverSettings();
+              
+                /*!
+                 * Retrieves whether a min/max equation solving technique has been set.
+                 *
+                 * @return True iff an equation solving technique has been set.
+                 */
+                bool isMinMaxEquationSolvingMethodSet() const;
+                
+                /*!
+                 * Retrieves the selected min/max equation solving technique.
+                 *
+                 * @return The selected min/max equation solving technique.
+                 */
+                storm::solver::MinMaxMethod getMinMaxEquationSolvingMethod() const;
+                
+                /*!
+                 * Retrieves whether the maximal iteration count has been set.
+                 *
+                 * @return True iff the maximal iteration count has been set.
+                 */
+                bool isMaximalIterationCountSet() const;
+                
+                /*!
+                 * Retrieves the maximal number of iterations to perform until giving up on converging.
+                 *
+                 * @return The maximal iteration count.
+                 */
+                uint_fast64_t getMaximalIterationCount() const;
+                
+                /*!
+                 * Retrieves whether the precision has been set.
+                 *
+                 * @return True iff the precision has been set.
+                 */
+                bool isPrecisionSet() const;
+                
+                /*!
+                 * Retrieves the precision that is used for detecting convergence.
+                 *
+                 * @return The precision to use for detecting convergence.
+                 */
+                double getPrecision() const;
+                
+                /*!
+                 * Retrieves whether the convergence criterion has been set.
+                 *
+                 * @return True iff the convergence criterion has been set.
+                 */
+                bool isConvergenceCriterionSet() const;
+                
+                /*!
+                 * Retrieves the selected convergence criterion.
+                 *
+                 * @return The selected convergence criterion.
+                 */
+                ConvergenceCriterion getConvergenceCriterion() const;
+                
+                // The name of the module.
+                static const std::string moduleName;
+                
+            private:
+                static const std::string solvingMethodOptionName;
+                static const std::string maximalIterationsOptionName;
+                static const std::string maximalIterationsOptionShortName;
+                static const std::string precisionOptionName;
+                static const std::string absoluteOptionName;
+            };
+            
+        }
+    }
+}

src/solver/EigenLinearEquationSolver.cpp

@@ -196,6 +196,34 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void EigenLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
+            // Typedef the map-type so we don't have to spell it out.
+            typedef decltype(Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
+
+            auto eigenX = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(x.data(), x.size());
+            auto eigenResult = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(result.data(), result.size());
+
+            std::unique_ptr<MapType> eigenB;
+            if (b != nullptr) {
+                eigenB = std::make_unique<MapType>(Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b->data(), b->size()));
+            }
+            
+            if (&x != &result) {
+                if (b != nullptr) {
+                    eigenResult.noalias() = *eigenA * eigenX + *eigenB;
+                } else {
+                    eigenResult.noalias() = *eigenA * eigenX;
+                }
+            } else {
+                if (b != nullptr) {
+                    eigenResult = *eigenA * eigenX + *eigenB;
+                } else {
+                    eigenResult = *eigenA * eigenX;
+                }
+            }
+        }
+        
         template<typename ValueType>
         void EigenLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
             // Typedef the map-type so we don't have to spell it out.
@@ -223,7 +251,6 @@ namespace storm {
                 } else {
                     nextX->noalias() = *eigenA * *currentX;
                 }
-                std::swap(nextX, currentX);
             }
             
             // If the last result we obtained is not the one in the input vector x, we swap the result there.
@@ -286,6 +313,7 @@ namespace storm {
                 } else {
                     nextX->noalias() = *eigenA * *currentX;
                 }
+                std::swap(nextX, currentX);
             }
             
             // If the last result we obtained is not the one in the input vector x, we swap the result there.
@@ -338,6 +366,7 @@ namespace storm {
                 } else {
                     nextX->noalias() = *eigenA * *currentX;
                 }
+                std::swap(nextX, currentX);
             }
             
             // If the last result we obtained is not the one in the input vector x, we swap the result there.

src/solver/EigenLinearEquationSolver.h

@@ -62,8 +62,8 @@ namespace storm {
             EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EigenLinearEquationSolverSettings<ValueType> const& settings = EigenLinearEquationSolverSettings<ValueType>());
             
             virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
 
             EigenLinearEquationSolverSettings<ValueType>& getSettings();
             EigenLinearEquationSolverSettings<ValueType> const& getSettings() const;

src/solver/EliminationLinearEquationSolver.cpp

@@ -137,6 +137,23 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void EliminationLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
+            if (&x != &result) {
+                A.multiplyWithVector(x, result);
+                if (b != nullptr) {
+                    storm::utility::vector::addVectors(result, *b, result);
+                }
+            } else {
+                // If the two vectors are aliases, we need to create a temporary.
+                std::vector<ValueType> tmp(result.size());
+                A.multiplyWithVector(x, tmp);
+                if (b != nullptr) {
+                    storm::utility::vector::addVectors(tmp, *b, result);
+                }
+            }
+        }
+        
         template<typename ValueType>
         EliminationLinearEquationSolverSettings<ValueType>& EliminationLinearEquationSolver<ValueType>::getSettings() {
             return settings;

src/solver/EliminationLinearEquationSolver.h

@@ -30,8 +30,8 @@ namespace storm {
             EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EliminationLinearEquationSolverSettings<ValueType> const& settings = EliminationLinearEquationSolverSettings<ValueType>());
             
             virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
 
             EliminationLinearEquationSolverSettings<ValueType>& getSettings();
             EliminationLinearEquationSolverSettings<ValueType> const& getSettings() const;

src/solver/GmmxxLinearEquationSolver.cpp

@@ -214,6 +214,14 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void GmmxxLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
+            gmm::mult(*gmmxxMatrix, x, result);
+            if (b != nullptr) {
+                gmm::add(*b, result);
+            }
+        }
+        
         template<typename ValueType>
         uint_fast64_t GmmxxLinearEquationSolver<ValueType>::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             // Get a Jacobi decomposition of the matrix A.

src/solver/GmmxxLinearEquationSolver.h

@@ -90,8 +90,8 @@ namespace storm {
             GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings = GmmxxLinearEquationSolverSettings<ValueType>());
             
             virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
 
             GmmxxLinearEquationSolverSettings<ValueType>& getSettings();
             GmmxxLinearEquationSolverSettings<ValueType> const& getSettings() const;

src/solver/LinearEquationSolver.h

@@ -40,13 +40,25 @@ namespace storm {
              * matrix A has to be given upon construction time of the solver object.
              *
              * @param x The initial vector with which to perform matrix-vector multiplication. Its length must be equal
-             * to the number of rows of A.
+             * to the number of columns of A.
              * @param b If non-null, this vector is added after each multiplication. If given, its length must be equal
              * to the number of rows of A.
              * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
              * vector must be equal to the number of rows of A.
              */
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const = 0;
+            
+            /*!
+             * Performs on matrix-vector multiplication x' = A*x + b.
+             *
+             * @param x The input vector with which to multiply the matrix. Its length must be equal
+             * to the number of columns of A.
+             * @param result The target vector into which to write the multiplication result. Its length must be equal
+             * to the number of rows 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.
+             */
+            virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const = 0;
         };
         
         template<typename ValueType>

src/solver/MinMaxLinearEquationSolver.cpp

@@ -7,7 +7,7 @@
 #include "src/solver/TopologicalMinMaxLinearEquationSolver.h"
 
 #include "src/settings/SettingsManager.h"
-#include "src/settings/modules/CoreSettings.h"
+#include "src/settings/modules/MinMaxEquationSolverSettings.h"
 
 #include "src/utility/macros.h"
 #include "src/exceptions/NotImplementedException.h"
@@ -118,10 +118,10 @@ namespace storm {
         template<typename MatrixType>
         std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> GeneralMinMaxLinearEquationSolverFactory<ValueType>::selectSolver(MatrixType&& matrix) const {
             std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> result;
-            auto technique = storm::settings::getModule<storm::settings::modules::CoreSettings>().getMinMaxEquationSolvingTechnique();
-            if (technique == MinMaxTechnique::ValueIteration || technique == MinMaxTechnique::PolicyIteration) {
+            auto method = storm::settings::getModule<storm::settings::modules::MinMaxEquationSolverSettings>().getMinMaxEquationSolvingMethod();
+            if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration) {
                 result = std::make_unique<StandardMinMaxLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix), std::make_unique<GeneralLinearEquationSolverFactory<ValueType>>());
-            } else if (technique == MinMaxTechnique::Topological) {
+            } else if (method == MinMaxMethod::Topological) {
                 result = std::make_unique<TopologicalMinMaxLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
             } else {
                 STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique.");

src/solver/NativeLinearEquationSolver.cpp

@@ -219,6 +219,23 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void NativeLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
+            if (&x != &result) {
+                A.multiplyWithVector(x, result);
+                if (b != nullptr) {
+                    storm::utility::vector::addVectors(result, *b, result);
+                }
+            } else {
+                // If the two vectors are aliases, we need to create a temporary.
+                std::vector<ValueType> tmp(result.size());
+                A.multiplyWithVector(x, tmp);
+                if (b != nullptr) {
+                    storm::utility::vector::addVectors(tmp, *b, result);
+                }
+            }
+        }
+        
         template<typename ValueType>
         NativeLinearEquationSolverSettings<ValueType>& NativeLinearEquationSolver<ValueType>::getSettings() {
             return settings;

src/solver/NativeLinearEquationSolver.h

@@ -47,8 +47,8 @@ namespace storm {
             NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, NativeLinearEquationSolverSettings<ValueType> const& settings = NativeLinearEquationSolverSettings<ValueType>());
             
             virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
 
             NativeLinearEquationSolverSettings<ValueType>& getSettings();
             NativeLinearEquationSolverSettings<ValueType> const& getSettings() const;

src/solver/SolverSelectionOptions.cpp

@@ -2,13 +2,13 @@
 
 namespace storm {
     namespace solver {
-        std::string toString(MinMaxTechnique m) {
+        std::string toString(MinMaxMethod m) {
             switch(m) {
-                case MinMaxTechnique::PolicyIteration:
+                case MinMaxMethod::PolicyIteration:
                     return "policy";
-                case MinMaxTechnique::ValueIteration:
+                case MinMaxMethod::ValueIteration:
                     return "value";
-                case MinMaxTechnique::Topological:
+                case MinMaxMethod::Topological:
                     return "topological";
 
             }

src/solver/SolverSelectionOptions.h

@@ -6,7 +6,7 @@
 
 namespace storm {
     namespace solver {
-        ExtendEnumsWithSelectionField(MinMaxTechnique, PolicyIteration, ValueIteration, Topological)
+        ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, Topological)
         
         ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk)
         ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination)

src/solver/StandardMinMaxLinearEquationSolver.cpp

@@ -1,13 +1,14 @@
 #include "src/solver/StandardMinMaxLinearEquationSolver.h"
 
 #include "src/settings/SettingsManager.h"
-#include "src/settings/modules/CoreSettings.h"
+#include "src/settings/modules/MinMaxEquationSolverSettings.h"
 
 #include "src/solver/GmmxxLinearEquationSolver.h"
 #include "src/solver/EigenLinearEquationSolver.h"
 #include "src/solver/NativeLinearEquationSolver.h"
 #include "src/solver/EliminationLinearEquationSolver.h"
 
+#include "src/utility/vector.h"
 #include "src/utility/macros.h"
 #include "src/exceptions/InvalidSettingsException.h"
 
@@ -15,10 +16,15 @@ namespace storm {
     namespace solver {
         
         StandardMinMaxLinearEquationSolverSettings::StandardMinMaxLinearEquationSolverSettings() {
-            auto technique = storm::settings::getModule<storm::settings::modules::CoreSettings>().getMinMaxEquationSolvingTechnique();
-            switch (technique) {
-                case MinMaxTechnique::ValueIteration: this->solutionMethod = SolutionMethod::ValueIteration; break;
-                case MinMaxTechnique::PolicyIteration: this->solutionMethod = SolutionMethod::PolicyIteration; break;
+            // Get the settings object to customize linear solving.
+            storm::settings::modules::MinMaxEquationSolverSettings const& settings = storm::settings::getModule<storm::settings::modules::MinMaxEquationSolverSettings>();
+
+            maximalNumberOfIterations = settings.getMaximalIterationCount();
+            
+            auto method = settings.getMinMaxEquationSolvingMethod();
+            switch (method) {
+                case MinMaxMethod::ValueIteration: this->solutionMethod = SolutionMethod::ValueIteration; break;
+                case MinMaxMethod::PolicyIteration: this->solutionMethod = SolutionMethod::PolicyIteration; break;
                 default:
                     STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique.");
             }
@@ -28,10 +34,34 @@ namespace storm {
             this->solutionMethod = solutionMethod;
         }
         
+        void StandardMinMaxLinearEquationSolverSettings::setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations) {
+            this->maximalNumberOfIterations = maximalNumberOfIterations;
+        }
+        
+        void StandardMinMaxLinearEquationSolverSettings::setRelativeTerminationCriterion(bool value) {
+            this->relative = value;
+        }
+        
+        void StandardMinMaxLinearEquationSolverSettings::setPrecision(double precision) {
+            this->precision = precision;
+        }
+        
         StandardMinMaxLinearEquationSolverSettings::SolutionMethod const& StandardMinMaxLinearEquationSolverSettings::getSolutionMethod() const {
             return solutionMethod;
         }
         
+        uint64_t StandardMinMaxLinearEquationSolverSettings::getMaximalNumberOfIterations() const {
+            return maximalNumberOfIterations;
+        }
+        
+        double StandardMinMaxLinearEquationSolverSettings::getPrecision() const {
+            return precision;
+        }
+        
+        bool StandardMinMaxLinearEquationSolverSettings::getRelativeTerminationCriterion() const {
+            return relative;
+        }
+        
         template<typename ValueType>
         StandardMinMaxLinearEquationSolver<ValueType>::StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings) : settings(settings), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), localA(nullptr), A(A) {
             // Intentionally left empty.
@@ -43,13 +73,98 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
-            //FIXME: implement
+        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+            switch (this->getSettings().getSolutionMethod()) {
+                case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration: solveEquationSystemValueIteration(dir, x, b, multiplyResult, newX); break;
+                case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: solveEquationSystemPolicyIteration(dir, x, b, multiplyResult, newX); break;
+            }
+        }
+        
+        template<typename ValueType>
+        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationSystemPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+            // FIXME.
+        }
+        
+        template<typename ValueType>
+        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationSystemValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+            // Create scratch memory if none was provided.
+            bool multiplyResultMemoryProvided = multiplyResult != nullptr;
+            if (multiplyResult == nullptr) {
+                multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
+            }
+            std::vector<ValueType>* currentX = &x;
+            bool xMemoryProvided = newX != nullptr;
+            if (newX == nullptr) {
+                newX = new std::vector<ValueType>(x.size());
+            }
+            
+            // Keep track of which of the vectors for x is the auxiliary copy.
+            std::vector<ValueType>* copyX = newX;
+            
+            uint64_t iterations = 0;
+            bool converged = false;
+            
+            // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations.
+            while (!converged && iterations < this->getSettings().getMaximalNumberOfIterations() && (!this->hasCustomTerminationCondition() || this->getTerminationCondition().terminateNow(*currentX))) {
+                // Compute x' = A*x + b.
+                this->A.multiplyWithVector(*currentX, *multiplyResult);
+                storm::utility::vector::addVectors(*multiplyResult, b, *multiplyResult);
+                
+                // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
+                // element of the min/max operator stack.
+                storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, *newX, this->A.getRowGroupIndices());
+                
+                // Determine whether the method converged.
+                converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion());
+                
+                // Update environment variables.
+                std::swap(currentX, newX);
+                ++iterations;
+            }
+            
+            // Check if the solver converged and issue a warning otherwise.
+            if (converged) {
+                STORM_LOG_INFO("Iterative solver converged after " << iterations << " iterations.");
+            } else {
+                STORM_LOG_WARN("Iterative solver did not converge after " << iterations << " iterations.");
+            }
+            
+            // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
+            // is currently stored in currentX, but x is the output vector.
+            if (currentX == copyX) {
+                std::swap(x, *currentX);
+            }
+            
+            if (!xMemoryProvided) {
+                delete copyX;
+            }
+            
+            if (!multiplyResultMemoryProvided) {
+                delete multiplyResult;
+            }
         }
         
         template<typename ValueType>
-        void StandardMinMaxLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
-            //FIXME: implement
+        void StandardMinMaxLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(A);
+            
+            // If scratch memory was not provided, we create it.
+            bool multiplyResultMemoryProvided = multiplyResult != nullptr;
+            if (!multiplyResult) {
+                multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
+            }
+            
+            for (uint64_t i = 0; i < n; ++i) {
+                solver->performMatrixVectorMultiplication(x, *multiplyResult, b);
+                
+                // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
+                // element of the min/max operator stack.
+                storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, this->A.getRowGroupIndices());
+            }
+            
+            if (!multiplyResultMemoryProvided) {
+                delete multiplyResult;
+            }
         }
         
         template<typename ValueType>

src/solver/StandardMinMaxLinearEquationSolver.h

@@ -15,11 +15,20 @@ namespace storm {
             };
             
             void setSolutionMethod(SolutionMethod const& solutionMethod);
+            void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations);
+            void setRelativeTerminationCriterion(bool value);
+            void setPrecision(double precision);
 
             SolutionMethod const& getSolutionMethod() const;
+            uint64_t getMaximalNumberOfIterations() const;
+            double getPrecision() const;
+            bool getRelativeTerminationCriterion() const;
 
         private:
             SolutionMethod solutionMethod;
+            uint64_t maximalNumberOfIterations;
+            double precision;
+            bool relative;
         };
         
         template<typename ValueType>
@@ -28,13 +37,16 @@ namespace storm {
             StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings());
             StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings());
             
-            virtual void solveEquationSystem(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
-            virtual void performMatrixVectorMultiplication(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            virtual void solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
 
             StandardMinMaxLinearEquationSolverSettings const& getSettings() const;
             StandardMinMaxLinearEquationSolverSettings& getSettings();
             
         private:
+            void solveEquationSystemPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const;
+            void solveEquationSystemValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const;
+
             /// The settings of this solver.
             StandardMinMaxLinearEquationSolverSettings settings;
             

test/functional/solver/MinMaxTechniqueSelectionTest.cpp

@@ -2,9 +2,9 @@
 
 #include "src/solver/MinMaxLinearEquationSolver.h"
 
-TEST( MinMaxTechnique, Simple ) {
-    storm::solver::MinMaxTechniqueSelection ts = storm::solver::MinMaxTechniqueSelection::PolicyIteration;
-    storm::solver::MinMaxTechnique t = storm::solver::MinMaxTechnique::PolicyIteration;
+TEST( MinMaxMethod, Simple ) {
+    storm::solver::MinMaxMethodSelection ts = storm::solver::MinMaxMethodSelection::PolicyIteration;
+    storm::solver::MinMaxMethod t = storm::solver::MinMaxMethod::PolicyIteration;
        ASSERT_EQ(convert(ts), t);