Started refactoring the linear equation system solvers.

Former-commit-id: 72d647fd42

src/adapters/GmmxxAdapter.h

@@ -31,12 +31,12 @@ public:
 	 * @return A pointer to a row-major sparse matrix in gmm++ format.
 	 */
 	template<class T>
-	static gmm::csr_matrix<T>* toGmmxxSparseMatrix(storm::storage::SparseMatrix<T> const& matrix) {
+	static std::unique_ptr<gmm::csr_matrix<T>> toGmmxxSparseMatrix(storm::storage::SparseMatrix<T> const& matrix) {
 		uint_fast64_t realNonZeros = matrix.getEntryCount();
 		LOG4CPLUS_DEBUG(logger, "Converting matrix with " << realNonZeros << " non-zeros to gmm++ format.");
 
 		// Prepare the resulting matrix.
-		gmm::csr_matrix<T>* result = new gmm::csr_matrix<T>(matrix.rowCount, matrix.columnCount);
+        std::unique_ptr<gmm::csr_matrix<T>> result(new gmm::csr_matrix<T>(matrix.getRowCount(), matrix.getColumnCount()));
 
 		// Copy Row Indications
 		std::copy(matrix.rowIndications.begin(), matrix.rowIndications.end(), result->jc.begin());
@@ -65,13 +65,13 @@ public:
 	 * @return A pointer to a row-major sparse matrix in gmm++ format.
 	 */
 	template<class T>
-	static gmm::csr_matrix<T>* toGmmxxSparseMatrix(storm::storage::SparseMatrix<T>&& matrix) {
+	static std::unique_ptr<gmm::csr_matrix<T>> toGmmxxSparseMatrix(storm::storage::SparseMatrix<T>&& matrix) {
 		uint_fast64_t realNonZeros = matrix.getEntryCount();
         std::cout << "here?!" << std::endl;
 		LOG4CPLUS_DEBUG(logger, "Converting matrix with " << realNonZeros << " non-zeros to gmm++ format.");
 
 		// Prepare the resulting matrix.
-		gmm::csr_matrix<T>* result = new gmm::csr_matrix<T>(matrix.rowCount, matrix.columnCount);
+        std::unique_ptr<gmm::csr_matrix<T>> result(new gmm::csr_matrix<T>(matrix.getRowCount(), matrix.getColumnCount()));
 
 		typedef unsigned long long IND_TYPE;
         typedef std::vector<IND_TYPE> vectorType_ull_t;

src/solver/AbstractLinearEquationSolver.h

@@ -1,23 +1,52 @@
 #ifndef STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_
 #define STORM_SOLVER_ABSTRACTLINEAREQUATIONSOLVER_H_
 
-#include "src/storage/SparseMatrix.h"
-
 #include <vector>
 
+#include "src/storage/SparseMatrix.h"
+
 namespace storm {
     namespace solver {
         
+        /*!
+         * A class that represents an abstract linear equation solver. In addition to solving a system of linear
+         * equations, the functionality to repeatedly multiply a matrix with a given vector is provided.
+         */
         template<class Type>
         class AbstractLinearEquationSolver {
         public:
-            
+            /*!
+             * Makes a copy of the linear equation solver.
+             *
+             * @return A pointer to a copy of the linear equation solver.
+             */
             virtual AbstractLinearEquationSolver<Type>* clone() const = 0;
             
-            virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const = 0;
-            
-            virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const = 0;
+            /*!
+             * Solves the equation system A*x = b. The matrix A is required to be square and have a unique solution.
+             * The solution of the set of linear equations will be written to the vector x.
+             *
+             * @param A The coefficient matrix of the equation system.
+             * @param x The solution vector that has to be computed. Its length must be equal to the number of rows of A.
+             * @param b The right-hand side of the equation system. 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 solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<Type>* multiplyResult = nullptr) const = 0;
             
+            /*!
+             * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After
+             * performing the necessary multiplications, the result is written to the input vector x.
+             *
+             * @param A The matrix to use for matrix-vector multiplication.
+             * @param x The initial vector with which to perform matrix-vector multiplication. Its length must be equal
+             * to the number of rows 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(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1, std::vector<Type>* multiplyResult = nullptr) const = 0;
         };
         
     } // namespace solver

src/solver/AbstractNondeterministicLinearEquationSolver.h

@@ -16,8 +16,8 @@ namespace storm {
             AbstractNondeterministicLinearEquationSolver() {
                 storm::settings::Settings* s = storm::settings::Settings::getInstance();
                 precision = s->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
-                maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger();
-                relative = s->getOptionByLongName("relative").getArgument(0).getValueAsBoolean();
+                maxIterations = s->getOptionByLongName("maxiter").getArgument(0).getValueAsUnsignedInteger();
+                relative = !s->isSet("absolute");
             }
             
             AbstractNondeterministicLinearEquationSolver(double precision, uint_fast64_t maxIterations, bool relative) : precision(precision), maxIterations(maxIterations), relative(relative) {

src/solver/GmmxxLinearEquationSolver.cpp

@@ -2,27 +2,277 @@
 
 #include <vector>
 #include <string>
+#include <cmath>
 
-bool GmmxxLinearEquationSolverOptionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* instance) -> bool {
+#include "src/adapters/GmmxxAdapter.h"
+#include "src/settings/Settings.h"
+#include "src/utility/vector.h"
+#include "src/utility/constants.h"
+#include "src/exceptions/InvalidStateException.h"
+
+#include "gmm/gmm_matrix.h"
+#include "gmm/gmm_iter_solvers.h"
 
+bool GmmxxLinearEquationSolverOptionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* instance) -> bool {
+    // Offer all available methods as a command line option.
 	std::vector<std::string> methods;
 	methods.push_back("bicgstab");
 	methods.push_back("qmr");
-	methods.push_back("lscg");
 	methods.push_back("gmres");
 	methods.push_back("jacobi");
+	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "gmmlin", "", "The method to be used for solving linear equation systems with the gmm++ engine. Available are: bicgstab, qmr, gmres, jacobi.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("gmres").build()).build());
 
-	std::vector<std::string> preconditioner;
+    // Register available preconditioners.
+    std::vector<std::string> preconditioner;
 	preconditioner.push_back("ilu");
 	preconditioner.push_back("diagonal");
 	preconditioner.push_back("ildlt");
 	preconditioner.push_back("none");
+	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "gmmpre", "", "The preconditioning technique used for solving linear equation systems with the gmm++ engine. Available are: ilu, diagonal, none.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the preconditioning method.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(preconditioner)).setDefaultValueString("ilu").build()).build());
 
-	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "leMethod", "", "The Method used in Linear Equation Solving. Available are: bicgstab, qmr, lscg, gmres, jacobi").addArgument(storm::settings::ArgumentBuilder::createStringArgument("leMethodName", "The Name of the Method to use").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("gmres").build()).build());
-	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "preconditioner", "", "The Preconditioning Technique used in Linear Equation Solving. Available are: ilu, diagonal, ildlt, none").addArgument(storm::settings::ArgumentBuilder::createStringArgument("preconditionerName", "The Name of the Preconditioning Method").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(preconditioner)).setDefaultValueString("ilu").build()).build());
-	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "maxIterations", "", "Maximum number of Iterations to perform while solving a linear equation system").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("iterationCount", "Max. Iteration Count").setDefaultValueUnsignedInteger(10000).build()).build());
-	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "precision", "", "Precision used for iterative solving of linear equation systems").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("precisionValue", "Precision").setDefaultValueDouble(1e-6).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1000000.0)).build()).build());
-	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "relative", "", "Whether the relative or the absolute error is considered for deciding convergence via a given precision").addArgument(storm::settings::ArgumentBuilder::createBooleanArgument("useRelative", "relative or absolute comparison").setDefaultValueBoolean(true).build()).build());
+    instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "gmmrestart", "", "The number of iteration until restarted methods are actually restarted.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The number of iterations.").setDefaultValueUnsignedInteger(50).build()).build());
+
+	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "maxiter", "i", "The maximal number of iterations to perform before iterative solving is aborted.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(10000).build()).build());
+    
+	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "precision", "", "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());
+    
+	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "absolute", "", "Whether the relative or the absolute error is considered for deciding convergence.").build());
 
 	return true;
-});
+});
+
+namespace storm {
+    namespace solver {
+        
+        template<typename ValueType>
+        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(SolutionMethod method, double precision, bool relative, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, uint_fast64_t restart) : method(method), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner), restart(restart) {
+            // Intentionally left empty.
+        }
+        
+
+        template<typename ValueType>
+        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver() {
+            // Get the settings object to customize linear solving.
+            storm::settings::Settings* settings = storm::settings::Settings::getInstance();
+            
+            // Get appropriate settings.
+            maximalNumberOfIterations = settings->getOptionByLongName("maxiter").getArgument(0).getValueAsUnsignedInteger();
+            precision = settings->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
+            relative = settings->isSet("absolute");
+            restart = settings->getOptionByLongName("gmmrestart").getArgument(0).getValueAsUnsignedInteger();
+            
+            // Determine the method to be used.
+            std::string const& methodAsString = settings->getOptionByLongName("gmmlin").getArgument(0).getValueAsString();
+            if (methodAsString == "bicgstab") {
+                method = BICGSTAB;
+            } else if (methodAsString == "qmr") {
+                method = QMR;
+            } else if (methodAsString == "gmres") {
+                method = GMRES;
+            } else if (methodAsString == "jacobi") {
+                method = JACOBI;
+            }
+            
+            // Check which preconditioner to use.
+            std::string const& preconditionAsString = settings->getOptionByLongName("gmmpre").getArgument(0).getValueAsString();
+            if (preconditionAsString == "ilu") {
+                preconditioner = ILU;
+            } else if (preconditionAsString == "diagonal") {
+                preconditioner = DIAGONAL;
+            } else if (preconditionAsString == "none") {
+                preconditioner = NONE;
+            }
+        }
+        
+        template<typename ValueType>
+        AbstractLinearEquationSolver<ValueType>* GmmxxLinearEquationSolver<ValueType>::clone() const {
+            return new GmmxxLinearEquationSolver<ValueType>(*this);
+        }
+        
+        template<typename ValueType>
+        void GmmxxLinearEquationSolver<ValueType>::solveEquationSystem(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+            LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner " << preconditionerToString() << "'.");
+            if (method == JACOBI && preconditioner != NONE) {
+                LOG4CPLUS_WARN(logger, "Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored.");
+            }
+            
+            if (method == BICGSTAB || method == QMR || method == GMRES) {
+                std::unique_ptr<gmm::csr_matrix<ValueType>> gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
+                
+                // Prepare an iteration object that determines the accuracy and the maximum number of iterations.
+                gmm::iteration iter(precision, 0, maximalNumberOfIterations);
+                
+                if (method == BICGSTAB) {
+                    if (preconditioner == ILU) {
+                        gmm::bicgstab(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmA), iter);
+                    } else if (preconditioner == DIAGONAL) {
+                        gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmA), iter);
+                    } else if (preconditioner == NONE) {
+                        gmm::bicgstab(*gmmA, x, b, gmm::identity_matrix(), iter);
+                    }
+                } else if (method == QMR) {
+                    if (preconditioner == ILU) {
+                        gmm::qmr(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmA), iter);
+                    } else if (preconditioner == DIAGONAL) {
+                        gmm::qmr(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmA), iter);
+                    } else if (preconditioner == NONE) {
+                        gmm::qmr(*gmmA, x, b, gmm::identity_matrix(), iter);
+                    }
+                } else if (method == GMRES) {
+                    if (preconditioner == ILU) {
+                        gmm::gmres(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmA), restart, iter);
+                    } else if (preconditioner == DIAGONAL) {
+                        gmm::gmres(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmA), restart, iter);
+                    } else if (preconditioner == NONE) {
+                        gmm::gmres(*gmmA, x, b, gmm::identity_matrix(), restart, iter);
+                    }
+                }
+                
+                // Check if the solver converged and issue a warning otherwise.
+                if (iter.converged()) {
+                    LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
+                } else {
+                    LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
+                }
+            } else if (method == JACOBI) {
+                uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b, multiplyResult);
+                
+                // Check if the solver converged and issue a warning otherwise.
+                if (iterations < maximalNumberOfIterations) {
+                    LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
+                } else {
+                    LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
+                }
+            }
+        }
+        
+        template<typename ValueType>
+        void GmmxxLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+            // Transform the transition probability A to the gmm++ format to use its arithmetic.
+            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
+            
+            // Set up some temporary variables so that we can just swap pointers instead of copying the result after
+            // each iteration.
+            std::vector<ValueType>* swap = nullptr;
+            std::vector<ValueType>* currentX = &x;
+            
+            bool multiplyResultProvided = true;
+            std::vector<ValueType>* nextX = multiplyResult;
+            if (nextX == nullptr) {
+                nextX = new std::vector<ValueType>(x.size());
+                multiplyResultProvided = false;
+            }
+            std::vector<ValueType> const* copyX = nextX;
+            
+            // Now perform matrix-vector multiplication as long as we meet the bound.
+            for (uint_fast64_t i = 0; i < n; ++i) {
+                gmm::mult(*gmmxxMatrix, *currentX, *nextX);
+                std::swap(nextX, currentX);
+                
+                // If requested, add an offset to the current result vector.
+                if (b != nullptr) {
+                    gmm::add(*b, *currentX);
+                }
+            }
+            
+            // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x,
+            // because the output is supposed to be stored in the input vector x.
+            if (currentX == copyX) {
+                std::swap(x, *currentX);
+            }
+            
+            // If the vector for the temporary multiplication result was not provided, we need to delete it.
+            if (multiplyResultProvided) {
+                delete copyX;
+            }
+        }
+        
+        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.
+            std::pair<storm::storage::SparseMatrix<ValueType>, storm::storage::SparseMatrix<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+            
+            // Convert the (inverted) diagonal matrix to gmm++'s format.
+            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmDinv = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.second));
+            // Convert the LU matrix to gmm++'s format.
+            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.first));
+            
+            // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
+            bool multiplyResultProvided = true;
+            std::vector<ValueType>* nextX = multiplyResult;
+            if (nextX == nullptr) {
+                nextX = new std::vector<ValueType>(x.size());
+                multiplyResultProvided = false;
+            }
+            std::vector<ValueType> const* copyX = nextX;
+            std::vector<ValueType>* currentX = &x;
+            
+            // Target vector for precision calculation.
+            std::vector<ValueType> tmpX(x.size());
+            
+            // Set up additional environment variables.
+            uint_fast64_t iterationCount = 0;
+            bool converged = false;
+            
+            while (!converged && iterationCount < maximalNumberOfIterations) {
+                // Compute D^-1 * (b - LU * x) and store result in nextX.
+                gmm::mult(*gmmLU, *currentX, tmpX);
+                gmm::add(b, gmm::scaled(tmpX, -storm::utility::constantOne<ValueType>()), tmpX);
+                gmm::mult(*gmmDinv, tmpX, *nextX);
+                
+                // Swap the two pointers as a preparation for the next iteration.
+                std::swap(nextX, currentX);
+                
+                // Now check if the process already converged within our precision.
+                converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, precision, relative);
+                
+                // Increase iteration count so we can abort if convergence is too slow.
+                ++iterationCount;
+            }
+            
+            // If the last iteration did not write to the original x we have to swap the contents, because the
+            // output has to be written to the input parameter x.
+            if (currentX == copyX) {
+                std::swap(x, *currentX);
+            }
+            
+            // If the vector for the temporary multiplication result was not provided, we need to delete it.
+            if (!multiplyResultProvided) {
+                delete copyX;
+            }
+            
+            return iterationCount;
+        }
+        
+        template<typename ValueType>
+        std::string GmmxxLinearEquationSolver<ValueType>::methodToString() const {
+            if (method == BICGSTAB) {
+                return "bicgstab";
+            } else if (method == QMR) {
+                return "qmr";
+            } else if (method == GMRES) {
+                return "gmres";
+            } else if (method == JACOBI) {
+                return "jacobi";
+            } else {
+                throw storm::exceptions::InvalidStateException() << "Illegal method '" << method << "' set in GmmxxLinearEquationSolver.";
+            }
+        }
+        
+        template<typename ValueType>
+        std::string GmmxxLinearEquationSolver<ValueType>::preconditionerToString() const {
+            if (preconditioner == ILU) {
+                return "ilu";
+            } else if (preconditioner == DIAGONAL) {
+                return "diagonal";
+            } else if (preconditioner == NONE) {
+                return "none";
+            } else {
+                throw storm::exceptions::InvalidStateException() << "Illegal preconditioner '" << preconditioner << "' set in GmmxxLinearEquationSolver.";
+            }
+        }
+        
+        // Explicitly instantiate the solver.
+        template class GmmxxLinearEquationSolver<double>;
+    }
+}

src/solver/GmmxxLinearEquationSolver.h

@@ -2,270 +2,100 @@
 #define STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_
 
 #include "AbstractLinearEquationSolver.h"
-#include "src/adapters/GmmxxAdapter.h"
-#include "src/utility/constants.h"
-#include "src/settings/Settings.h"
-#include "src/utility/vector.h"
-
-#include "gmm/gmm_matrix.h"
-#include "gmm/gmm_iter_solvers.h"
-
-#include <cmath>
 
 namespace storm {
     namespace solver {
 
-        template<class Type>
-        class GmmxxLinearEquationSolver : public AbstractLinearEquationSolver<Type> {
+        /*!
+         * A class that uses the gmm++ library to implement the AbstractLinearEquationSolver interface.
+         */
+        template<typename ValueType>
+        class GmmxxLinearEquationSolver : public AbstractLinearEquationSolver<ValueType> {
         public:
+            // An enumeration specifying the available preconditioners.
+            enum Preconditioner {
+                ILU, DIAGONAL, NONE
+            };
             
-            virtual AbstractLinearEquationSolver<Type>* clone() const {
-                return new GmmxxLinearEquationSolver<Type>();
-            }
+            // An enumeration specifying the available solution methods.
+            enum SolutionMethod {
+                BICGSTAB, QMR, GMRES, JACOBI
+            };
             
-            virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const {
-                // Get the settings object to customize linear solving.
-                storm::settings::Settings* s = storm::settings::Settings::getInstance();
-                
-                // Prepare an iteration object that determines the accuracy, maximum number of iterations
-                // and the like.
-				uint_fast64_t maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger();
-                gmm::iteration iter(s->getOptionByLongName("precision").getArgument(0).getValueAsDouble(), 0, maxIterations);
-                
-                // Print some information about the used preconditioner.
-                std::string const precond = s->getOptionByLongName("preconditioner").getArgument(0).getValueAsString();
-                LOG4CPLUS_INFO(logger, "Starting iterative solver.");
-
-				// ALL available solvers must be declared in the cpp File, where the options are registered!
-				// Dito for the Preconditioners
-
-				std::string const chosenLeMethod = s->getOptionByLongName("leMethod").getArgument(0).getValueAsString();
-                if (chosenLeMethod == "jacobi") {
-                    if (precond != "none") {
-                        LOG4CPLUS_WARN(logger, "Requested preconditioner '" << precond << "', which is unavailable for the Jacobi method. Dropping preconditioner.");
-                    }
-                } else {
-                    if (precond == "ilu") {
-                        LOG4CPLUS_INFO(logger, "Using ILU preconditioner.");
-                    } else if (precond == "diagonal") {
-                        LOG4CPLUS_INFO(logger, "Using diagonal preconditioner.");
-                    } else if (precond == "ildlt") {
-                        LOG4CPLUS_INFO(logger, "Using ILDLT preconditioner.");
-                    } else if (precond == "none") {
-                        LOG4CPLUS_INFO(logger, "Using no preconditioner.");
-                    }
-                }
-                
-                // Now do the actual solving.
-                if (chosenLeMethod == "bicgstab") {
-                    LOG4CPLUS_INFO(logger, "Using BiCGStab method.");
-                    // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-                    gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
-                    if (precond == "ilu") {
-                        gmm::bicgstab(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
-                    } else if (precond == "diagonal") {
-                        gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
-                    } else if (precond == "ildlt") {
-                        gmm::bicgstab(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
-                    } else if (precond == "none") {
-                        gmm::bicgstab(*gmmA, x, b, gmm::identity_matrix(), iter);
-                    }
-                    
-                    // Check if the solver converged and issue a warning otherwise.
-                    if (iter.converged()) {
-                        LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
-                    } else {
-                        LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
-                    }
-                    delete gmmA;
-                } else if (chosenLeMethod == "qmr") {
-                    LOG4CPLUS_INFO(logger, "Using QMR method.");
-                    // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-                    gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
-                    if (precond == "ilu") {
-                        gmm::qmr(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
-                    } else if (precond == "diagonal") {
-                        gmm::qmr(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
-                    } else if (precond == "ildlt") {
-                        gmm::qmr(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*gmmA), iter);
-                    } else if (precond == "none") {
-                        gmm::qmr(*gmmA, x, b, gmm::identity_matrix(), iter);
-                    }
-                    
-                    // Check if the solver converged and issue a warning otherwise.
-                    if (iter.converged()) {
-                        LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
-                    } else {
-                        LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
-                    }
-                    delete gmmA;
-                } else if (chosenLeMethod == "lscg") {
-                    LOG4CPLUS_INFO(logger, "Using LSCG method.");
-                    // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-                    gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
-                    
-                    if (precond != "none") {
-                        LOG4CPLUS_WARN(logger, "Requested preconditioner '" << precond << "', which is unavailable for the LSCG method. Dropping preconditioner.");
-                    }
-                    gmm::least_squares_cg(*gmmA, x, b, iter);
-                    
-                    // Check if the solver converged and issue a warning otherwise.
-                    if (iter.converged()) {
-                        LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
-                    } else {
-                        LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
-                    }
-                    delete gmmA;
-                } else if (chosenLeMethod == "gmres") {
-                    LOG4CPLUS_INFO(logger, "Using GMRES method.");
-                    // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-                    gmm::csr_matrix<Type>* gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
-                    if (precond == "ilu") {
-                        gmm::gmres(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*gmmA), 50, iter);
-                    } else if (precond == "diagonal") {
-                        gmm::gmres(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*gmmA), 50, iter);
-                    } else if (precond == "ildlt") {
-                        gmm::gmres(*gmmA, x, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*gmmA), 50, iter);
-                    } else if (precond == "none") {
-                        gmm::gmres(*gmmA, x, b, gmm::identity_matrix(), 50, iter);
-                    }
-                    
-                    // Check if the solver converged and issue a warning otherwise.
-                    if (iter.converged()) {
-                        LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
-                    } else {
-                        LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
-                    }
-                    delete gmmA;
-                } else if (chosenLeMethod == "jacobi") {
-                    LOG4CPLUS_INFO(logger, "Using Jacobi method.");
-                    uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b);
-                    uint_fast64_t maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger();
-                    // Check if the solver converged and issue a warning otherwise.
-					if (iterations < maxIterations) {
-                        LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
-                    } else {
-                        LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
-                    }
-                }
-            }
+            /*!
+             * Constructs a linear equation solver with the given parameters.
+             *
+             * @param method The method to use for linear equation solving.
+             * @param precision The precision to use for convergence detection.
+             * @param relative If set, the relative error rather than the absolute error is considered for convergence
+             * detection.
+             * @param maximalNumberOfIterations The maximal number of iterations do perform before iteration is aborted.
+             * @param preconditioner The preconditioner to use when solving linear equation systems.
+             * @param restart An optional argument that specifies after how many iterations restarted methods are
+             * supposed to actually to a restart.
+             */
+            GmmxxLinearEquationSolver(SolutionMethod method, double precision, bool relative, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, uint_fast64_t restart = 0);
             
-            virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b, uint_fast64_t n = 1) const {
-                // Transform the transition probability A to the gmm++ format to use its arithmetic.
-                gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
-                
-                // Set up some temporary variables so that we can just swap pointers instead of copying the result after each
-                // iteration.
-                std::vector<Type>* swap = nullptr;
-                std::vector<Type>* currentVector = &x;
-                std::vector<Type>* tmpVector = new std::vector<Type>(A.getRowCount());
-                
-                // Now perform matrix-vector multiplication as long as we meet the bound.
-                for (uint_fast64_t i = 0; i < n; ++i) {
-                    gmm::mult(*gmmxxMatrix, *currentVector, *tmpVector);
-                    swap = tmpVector;
-                    tmpVector = currentVector;
-                    currentVector = swap;
-                    
-                    // If requested, add an offset to the current result vector.
-                    if (b != nullptr) {
-                        gmm::add(*b, *currentVector);
-                    }
-                }
-                
-                // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, because
-                // the output is supposed to be stored in x.
-                if (n % 2 == 1) {
-                    std::swap(x, *currentVector);
-                    delete currentVector;
-                } else {
-                    delete tmpVector;
-                }
-                
-                delete gmmxxMatrix;
-            }
+            /*!
+             * Constructs a linear equation solver with parameters being set according to the settings object.
+             */
+            GmmxxLinearEquationSolver();
+            
+            virtual AbstractLinearEquationSolver<ValueType>* clone() const override;
+            
+            virtual void solveEquationSystem(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            
+            virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
         private:
             /*!
              * Solves the linear equation system A*x = b given by the parameters using the Jacobi method.
              *
              * @param A The matrix specifying the coefficients of the linear equations.
-             * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but
-             * may be ignored.
+             * @param x The solution vector x. The initial values of x represent a guess of the real values to the
+             * solver, but may be set to zero.
              * @param b The right-hand side of the equation system.
-             * @returns The solution vector x of the system of linear equations as the content of the parameter x.
-             * @returns The number of iterations needed until convergence.
+             * @return The number of iterations needed until convergence if the solver converged and
+             * maximalNumberOfIteration otherwise.
+             * @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.
+             */
+            uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const;
+
+            /*!
+             * Retrieves the string representation of the solution method associated with this solver.
+             *
+             * @return The string representation of the solution method associated with this solver.
+             */
+            std::string methodToString() const;
+            
+            /*!
+             * Retrieves the string representation of the preconditioner associated with this solver.
+             *
+             * @return The string representation of the preconditioner associated with this solver.
              */
-            uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const {
-                // Get the settings object to customize linear solving.
-                storm::settings::Settings* s = storm::settings::Settings::getInstance();
-                
-                double precision = s->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
-				uint_fast64_t maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger();
-				bool relative = s->getOptionByLongName("relative").getArgument(0).getValueAsBoolean();
-                
-                // Get a Jacobi decomposition of the matrix A.
-                std::pair<storm::storage::SparseMatrix<Type>, storm::storage::SparseMatrix<Type>> jacobiDecomposition = A.getJacobiDecomposition();
-                
-                // Convert the (inverted) diagonal matrix to gmm++'s format.
-                gmm::csr_matrix<Type>* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(std::move(jacobiDecomposition.second));
-                // Convert the LU matrix to gmm++'s format.
-                gmm::csr_matrix<Type>* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(std::move(jacobiDecomposition.first));
-                
-                LOG4CPLUS_INFO(logger, "Starting iterative Jacobi Solver.");
-                
-                // x_(k + 1) = D^-1 * (b  - R * x_k)
-                // In x we keep a copy of the result for swapping in the loop (e.g. less copy-back).
-                std::vector<Type>* xNext = new std::vector<Type>(x.size());
-                const std::vector<Type>* xCopy = xNext;
-                std::vector<Type>* xCurrent = &x;
-                
-                // Target vector for precision calculation.
-                std::vector<Type> tmpX(x.size());
-                
-                // Set up additional environment variables.
-                uint_fast64_t iterationCount = 0;
-                bool converged = false;
-                
-                while (!converged && iterationCount < maxIterations) {
-                    // R * x_k (xCurrent is x_k) -> xNext
-                    gmm::mult(*gmmxxLU, *xCurrent, tmpX);
-                    // b - R * x_k (xNext contains R * x_k) -> xNext
-                    gmm::add(b, gmm::scaled(tmpX, -1.0), tmpX);
-                    // D^-1 * (b - R * x_k) -> xNext
-                    gmm::mult(*gmmxxDiagonalInverted, tmpX, *xNext);
-                    
-                    // Swap xNext with xCurrent so that the next iteration can use xCurrent again without having to copy the
-                    // vector.
-                    std::vector<Type> *const swap = xNext;
-                    xNext = xCurrent;
-                    xCurrent = swap;
-                    
-                    // Now check if the process already converged within our precision.
-                    converged = storm::utility::vector::equalModuloPrecision(*xCurrent, *xNext, precision, relative);
-                    
-                    // Increase iteration count so we can abort if convergence is too slow.
-                    ++iterationCount;
-                }
-                
-                // If the last iteration did not write to the original x we have to swap the contents, because the output has to
-                // be written to the parameter x.
-                if (xCurrent == xCopy) {
-                    std::swap(x, *xCurrent);
-                }
-                
-                // As xCopy always points to the copy of x used for swapping, we can safely delete it.
-                delete xCopy;
-                
-                // Also delete the other dynamically allocated variables.
-                delete gmmxxDiagonalInverted;
-                delete gmmxxLU;
-                
-                return iterationCount;
-            }
+            std::string preconditionerToString() const;
             
+            // The method to use for solving linear equation systems.
+            SolutionMethod method;
+
+            // The required precision for the iterative methods.
+            double precision;
+
+            // Sets whether the relative or absolute error is to be considered for convergence detection. Not that this
+            // only applies to the Jacobi method for this solver.
+            bool relative;
+
+            // The maximal number of iterations to do before iteration is aborted.
+            uint_fast64_t maximalNumberOfIterations;
+
+            // The preconditioner to use when solving the linear equation system.
+            Preconditioner preconditioner;
+
+            // A restart value that determines when restarted methods shall do so.
+            uint_fast64_t restart;
         };
-        
+
     } // namespace solver
 } // namespace storm
 

src/solver/GmmxxNondeterministicLinearEquationSolver.cpp

@@ -0,0 +1,11 @@
+#include "src/settings/Settings.h"
+
+bool GmmxxNondeterministicLinearEquationSolverOptionsRegistered = storm::settings::Settings::registerNewModule([] (storm::settings::Settings* instance) -> bool {
+	instance->addOption(storm::settings::OptionBuilder("GmmxxNondeterminsticLinearEquationSolver", "maxiter", "i", "The maximal number of iterations to perform before iterative solving is aborted.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(10000).build()).build());
+    
+	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "precision", "", "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-6).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+    
+	instance->addOption(storm::settings::OptionBuilder("GmmxxLinearEquationSolver", "absolute", "", "Whether the relative or the absolute error is considered for deciding convergence.").build());
+    
+	return true;
+});

src/solver/GmmxxNondeterministicLinearEquationSolver.h

@@ -22,7 +22,7 @@ namespace storm {
             
             virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const override {
                 // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-                gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
+                std::unique_ptr<gmm::csr_matrix<Type>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
                 
                 // Create vector for result of multiplication, which is reduced to the result vector after
                 // each multiplication.
@@ -42,14 +42,11 @@ namespace storm {
                         storm::utility::vector::reduceVectorMax(multiplyResult, x, nondeterministicChoiceIndices);
                     }
                 }
-                
-                // Delete intermediate results and return result.
-                delete gmmxxMatrix;
             }
             
             virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* multiplyResult = nullptr, std::vector<Type>* newX = nullptr) const override {
                 // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-                gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
+                std::unique_ptr<gmm::csr_matrix<Type>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
                 
                 // Set up the environment for the power method.
                 bool multiplyResultMemoryProvided = true;
@@ -102,7 +99,6 @@ namespace storm {
                 } else if (!xMemoryProvided) {
                     delete newX;
                 }
-                delete gmmxxMatrix;
         
                 if (!multiplyResultMemoryProvided) {
                     delete multiplyResult;