fixed some TBB-related issues and added power method for linear equation systems

src/storm/settings/modules/NativeEquationSolverSettings.cpp

@@ -24,7 +24,7 @@ namespace storm {
             const std::string NativeEquationSolverSettings::absoluteOptionName = "absolute";
             
             NativeEquationSolverSettings::NativeEquationSolverSettings() : ModuleSettings(moduleName) {
-                std::vector<std::string> methods = { "jacobi", "gaussseidel", "sor", "walkerchae" };
+                std::vector<std::string> methods = { "jacobi", "gaussseidel", "sor", "walkerchae", "power" };
                 this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the native engine.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(methods)).setDefaultValueString("jacobi").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());
@@ -54,6 +54,8 @@ namespace storm {
                     return NativeEquationSolverSettings::LinearEquationMethod::SOR;
                 } else if (linearEquationSystemTechniqueAsString == "walkerchae") {
                     return NativeEquationSolverSettings::LinearEquationMethod::WalkerChae;
+                } else if (linearEquationSystemTechniqueAsString == "power") {
+                    return NativeEquationSolverSettings::LinearEquationMethod::Power;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected.");
             }

src/storm/settings/modules/NativeEquationSolverSettings.h

@@ -13,7 +13,7 @@ namespace storm {
             class NativeEquationSolverSettings : public ModuleSettings {
             public:
                 // An enumeration of all available methods for solving linear equations.
-                enum class LinearEquationMethod { Jacobi, GaussSeidel, SOR, WalkerChae };
+                enum class LinearEquationMethod { Jacobi, GaussSeidel, SOR, WalkerChae, Power };
                 
                 // An enumeration of all available convergence criteria.
                 enum class ConvergenceCriterion { Absolute, Relative };

src/storm/solver/GmmxxMultiplier.cpp

@@ -136,6 +136,7 @@ namespace storm {
 #endif
         }
         
+#ifdef STORM_HAVE_INTELTBB
         template<typename T>
         class TbbMultAddReduceFunctor {
         public:
@@ -199,10 +200,16 @@ namespace storm {
             std::vector<T>& result;
             std::vector<uint64_t>* choices;
         };
+#endif
         
         template<typename T>
         void GmmxxMultiplier<T>::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, gmm::csr_matrix<T> const& matrix, std::vector<T> const& x, std::vector<T> const* b, std::vector<T>& result, std::vector<uint64_t>* choices) const {
+#ifdef STORM_HAVE_INTELTBB
             tbb::parallel_for(tbb::blocked_range<unsigned long>(0, rowGroupIndices.size() - 1, 10), TbbMultAddReduceFunctor<T>(dir, rowGroupIndices, matrix, x, b, result, choices));
+#else
+            STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version.");
+            multAddReduce(dir, rowGroupIndices, matrix, x, b, result, choices);
+#endif
         }
         
         template<>

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -26,6 +26,8 @@ namespace storm {
                 method = SolutionMethod::SOR;
             } else if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::WalkerChae) {
                 method = SolutionMethod::WalkerChae;
+            } else if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Power) {
+                method = SolutionMethod::Power;
             } else {
                 STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "The selected solution technique is invalid for this solver.");
             }
@@ -321,21 +323,77 @@ namespace storm {
                 std::swap(x, *currentX);
             }
 
+            // Resize the solution to the right size.
+            x.resize(this->A->getRowCount());
+            
+            // Finalize solution vector.
+            storm::utility::vector::applyPointwise(x, x, [this] (ValueType const& value) { return value - walkerChaeData->t; } );
+            
+            if (!this->isCachingEnabled()) {
+                clearCache();
+            }
+
             if (converged) {
                 STORM_LOG_INFO("Iterative solver converged in " << iterations << " iterations.");
             } else {
                 STORM_LOG_WARN("Iterative solver did not converge in " << iterations << " iterations.");
             }
+
+            return converged;
+        }
+        
+        template<typename ValueType>
+        bool NativeLinearEquationSolver<ValueType>::solveEquationsPower(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+            // FIXME: This solver will not work for all input systems. More concretely, the current implementation will
+            // not work for systems that have a 0 on the diagonal. This is not a restriction of this technique in general
+            // but arbitrary matrices require pivoting, which is not currently implemented.
+            STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (Power)");
             
-            // Resize the solution to the right size.
-            x.resize(this->A->getRowCount());
+            // We need to revert the transformation into an equation system matrix, because the elimination procedure
+            // and the distance computation is based on the probability matrix instead.
+            storm::storage::SparseMatrix<ValueType> locallyConvertedMatrix;
+            if (localA) {
+                localA->convertToEquationSystem();
+            } else {
+                locallyConvertedMatrix = *A;
+                locallyConvertedMatrix.convertToEquationSystem();
+            }
+            storm::storage::SparseMatrix<ValueType> const& transitionMatrix = localA ? *localA : locallyConvertedMatrix;
             
-            // Finalize solution vector.
-            storm::utility::vector::applyPointwise(x, x, [this] (ValueType const& value) { return value - walkerChaeData->t; } );
+            if (!this->cachedRowVector) {
+                this->cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount());
+            }
+            
+            std::vector<ValueType>* currentX = &x;
+            std::vector<ValueType>* nextX = this->cachedRowVector.get();
+
+            bool converged = false;
+            uint64_t iterations = 0;
+            while (!converged && iterations < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
+                this->multiplier.multAdd(transitionMatrix, *currentX, &b, *nextX);
+                
+                // Now check if the process already converged within our precision.
+                converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion());
+                
+                // Set up next iteration.
+                std::swap(currentX, nextX);
+                ++iterations;
+            }
+            
+            if (currentX == this->cachedRowVector.get()) {
+                std::swap(x, *nextX);
+            }
             
             if (!this->isCachingEnabled()) {
                 clearCache();
             }
+            
+            if (converged) {
+                STORM_LOG_INFO("Iterative solver converged in " << iterations << " iterations.");
+            } else {
+                STORM_LOG_WARN("Iterative solver did not converge in " << iterations << " iterations.");
+            }
+
             return converged;
         }
         
@@ -347,6 +405,8 @@ namespace storm {
                 return this->solveEquationsJacobi(x, b);
             } else if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::WalkerChae) {
                 return this->solveEquationsWalkerChae(x, b);
+            } else if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::Power) {
+                return this->solveEquationsPower(x, b);
             }
             
             STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unknown solving technique.");

src/storm/solver/NativeLinearEquationSolver.h

@@ -14,7 +14,7 @@ namespace storm {
         class NativeLinearEquationSolverSettings {
         public:
             enum class SolutionMethod {
-                Jacobi, GaussSeidel, SOR, WalkerChae
+                Jacobi, GaussSeidel, SOR, WalkerChae, Power
             };
 
             NativeLinearEquationSolverSettings();
@@ -70,6 +70,7 @@ namespace storm {
             virtual bool solveEquationsSOR(std::vector<ValueType>& x, std::vector<ValueType> const& b, ValueType const& omega) const;
             virtual bool solveEquationsJacobi(std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
             virtual bool solveEquationsWalkerChae(std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            virtual bool solveEquationsPower(std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
             
             // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
             // when the solver is destructed.

src/storm/solver/NativeMultiplier.cpp

@@ -85,7 +85,7 @@ namespace storm {
             matrix.multiplyAndReduceParallel(dir, rowGroupIndices, x, b, result, choices);
 #else
             STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version.");
-            multAddReduce(dir, rowGroupIndices, x, b, result, choices);
+            multAddReduce(dir, rowGroupIndices, matrix, x, b, result, choices);
 #endif
         }
 

src/storm/utility/VectorHelper.cpp

@@ -20,7 +20,7 @@ namespace storm {
         template<typename ValueType>
         VectorHelper<ValueType>::VectorHelper() : doParallelize(storm::settings::getModule<storm::settings::modules::CoreSettings>().isUseIntelTbbSet()) {
 #ifndef STORM_HAVE_INTELTBB
-            STORM_LOG_THROW(!parallelize, storm::exceptions::InvalidSettingsException, "Cannot parallelize without TBB.");
+            STORM_LOG_THROW(!doParallelize, storm::exceptions::InvalidSettingsException, "Cannot parallelize without TBB.");
 #endif
         }