added tests for Eigen solver

Former-commit-id: ede9efcee2

src/adapters/EigenAdapter.cpp

@@ -0,0 +1 @@
+#include "src/adapters/EigenAdapter.h"

src/adapters/EigenAdapter.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <memory>
+
+#include "src/utility/eigen.h"
+
+#include "src/storage/SparseMatrix.h"
+
+namespace storm {
+    namespace adapters {
+        
+        class EigenAdapter {
+        public:
+            /*!
+             * Converts a sparse matrix into a sparse matrix in the gmm++ format.
+             * @return A pointer to a row-major sparse matrix in gmm++ format.
+             */
+            template<class ValueType>
+            static std::unique_ptr<Eigen::SparseMatrix<ValueType>> toEigenSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix) {
+                std::vector<Eigen::Triplet<ValueType>> triplets;
+                triplets.reserve(matrix.getNonzeroEntryCount());
+                
+                for (uint64_t row = 0; row < matrix.getRowCount(); ++row) {
+                    for (auto const& element : matrix.getRow(row)) {
+                        triplets.emplace_back(row, element.getColumn(), element.getValue());
+                    }
+                }
+                
+                std::unique_ptr<Eigen::SparseMatrix<ValueType>> result = std::make_unique<Eigen::SparseMatrix<ValueType>>(matrix.getRowCount(), matrix.getColumnCount());
+                result->setFromTriplets(triplets.begin(), triplets.end());
+                return result;
+            }
+        };
+        
+    }
+}

src/adapters/GmmxxAdapter.h

@@ -1,68 +1,59 @@
-/*
- * GmmxxAdapter.h
- *
- *  Created on: 24.12.2012
- *      Author: Christian Dehnert
- */
-
 #ifndef STORM_ADAPTERS_GMMXXADAPTER_H_
 #define STORM_ADAPTERS_GMMXXADAPTER_H_
 
-#include <new>
-#include <utility>
+#include <algorithm>
+#include <memory>
 
 #include "src/utility/gmm.h"
 
 #include "src/storage/SparseMatrix.h"
-#include "src/utility/ConversionHelper.h"
 
 #include "src/utility/macros.h"
 
 namespace storm {
     
-namespace adapters {
-
-class GmmxxAdapter {
-public:
-	/*!
-	 * Converts a sparse matrix into a sparse matrix in the gmm++ format.
-	 * @return A pointer to a row-major sparse matrix in gmm++ format.
-	 */
-	template<class T>
-	static std::unique_ptr<gmm::csr_matrix<T>> toGmmxxSparseMatrix(storm::storage::SparseMatrix<T> const& matrix) {
-		uint_fast64_t realNonZeros = matrix.getEntryCount();
-		STORM_LOG_DEBUG("Converting matrix with " << realNonZeros << " non-zeros to gmm++ format.");
+    namespace adapters {
         
-		// Prepare the resulting matrix.
-        std::unique_ptr<gmm::csr_matrix<T>> result(new gmm::csr_matrix<T>(matrix.getRowCount(), matrix.getColumnCount()));
+        class GmmxxAdapter {
+        public:
+            /*!
+             * Converts a sparse matrix into a sparse matrix in the gmm++ format.
+             * @return A pointer to a row-major sparse matrix in gmm++ format.
+             */
+            template<class T>
+            static std::unique_ptr<gmm::csr_matrix<T>> toGmmxxSparseMatrix(storm::storage::SparseMatrix<T> const& matrix) {
+                uint_fast64_t realNonZeros = matrix.getEntryCount();
+                STORM_LOG_DEBUG("Converting matrix with " << realNonZeros << " non-zeros to gmm++ format.");
                 
-		// Copy Row Indications
-		std::copy(matrix.rowIndications.begin(), matrix.rowIndications.end(), result->jc.begin());
+                // Prepare the resulting matrix.
+                std::unique_ptr<gmm::csr_matrix<T>> result(new gmm::csr_matrix<T>(matrix.getRowCount(), matrix.getColumnCount()));
                 
-        // Copy columns and values.
-        std::vector<T> values;
-        values.reserve(matrix.getEntryCount());
+                // Copy Row Indications
+                std::copy(matrix.rowIndications.begin(), matrix.rowIndications.end(), result->jc.begin());
                 
-        // To match the correct vector type for gmm, we create the vector with the exact same type.
-        decltype(result->ir) columns;
-        columns.reserve(matrix.getEntryCount());
+                // Copy columns and values.
+                std::vector<T> values;
+                values.reserve(matrix.getEntryCount());
                 
-        for (auto const& entry : matrix) {
-            columns.emplace_back(entry.getColumn());
-            values.emplace_back(entry.getValue());
-        }
+                // To match the correct vector type for gmm, we create the vector with the exact same type.
+                decltype(result->ir) columns;
+                columns.reserve(matrix.getEntryCount());
                 
-        std::swap(result->ir, columns);
-        std::swap(result->pr, values);
+                for (auto const& entry : matrix) {
+                    columns.emplace_back(entry.getColumn());
+                    values.emplace_back(entry.getValue());
+                }
                 
-		STORM_LOG_DEBUG("Done converting matrix to gmm++ format.");
+                std::swap(result->ir, columns);
+                std::swap(result->pr, values);
                 
-		return result;
-	}
-};
+                STORM_LOG_DEBUG("Done converting matrix to gmm++ format.");
                 
-} //namespace adapters
+                return result;
+            }
+        };
         
-} //namespace storm
+    } // namespace adapters
+} // namespace storm
 
 #endif /* STORM_ADAPTERS_GMMXXADAPTER_H_ */

src/settings/SettingsManager.cpp

@@ -20,6 +20,7 @@
 #include "src/settings/modules/CounterexampleGeneratorSettings.h"
 #include "src/settings/modules/CuddSettings.h"
 #include "src/settings/modules/SylvanSettings.h"
+#include "src/settings/modules/EigenEquationSolverSettings.h"
 #include "src/settings/modules/GmmxxEquationSolverSettings.h"
 #include "src/settings/modules/NativeEquationSolverSettings.h"
 #include "src/settings/modules/BisimulationSettings.h"
@@ -508,6 +509,7 @@ namespace storm {
             storm::settings::addModule<storm::settings::modules::CuddSettings>();
             storm::settings::addModule<storm::settings::modules::SylvanSettings>();
             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::BisimulationSettings>();
             storm::settings::addModule<storm::settings::modules::GlpkSettings>();

src/settings/modules/EigenEquationSolverSettings.cpp

@@ -0,0 +1,93 @@
+#include "src/settings/modules/EigenEquationSolverSettings.h"
+
+#include "src/settings/Option.h"
+#include "src/settings/OptionBuilder.h"
+#include "src/settings/ArgumentBuilder.h"
+#include "src/settings/Argument.h"
+
+#include "src/settings/SettingsManager.h"
+#include "src/settings/modules/MarkovChainSettings.h"
+#include "src/solver/SolverSelectionOptions.h"
+
+namespace storm {
+    namespace settings {
+        namespace modules {
+            
+            const std::string EigenEquationSolverSettings::moduleName = "eigen";
+            const std::string EigenEquationSolverSettings::techniqueOptionName = "method";
+            const std::string EigenEquationSolverSettings::preconditionOptionName = "precond";
+            const std::string EigenEquationSolverSettings::maximalIterationsOptionName = "maxiter";
+            const std::string EigenEquationSolverSettings::maximalIterationsOptionShortName = "i";
+            const std::string EigenEquationSolverSettings::precisionOptionName = "precision";
+            
+            EigenEquationSolverSettings::EigenEquationSolverSettings() : ModuleSettings(moduleName) {
+                std::vector<std::string> methods = {"sparselu", "bicgstab"};
+                this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the eigen solver. Available are {sparselu, bicgstab}.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("sparselu").build()).build());
+                
+                // Register available preconditioners.
+                std::vector<std::string> preconditioner = {"ilu", "diagonal", "none"};
+                this->addOption(storm::settings::OptionBuilder(moduleName, preconditionOptionName, true, "The preconditioning technique used for solving linear equation systems. 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());
+                
+                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());
+            }
+            
+            bool EigenEquationSolverSettings::isLinearEquationSystemMethodSet() const {
+                return this->getOption(techniqueOptionName).getHasOptionBeenSet();
+            }
+            
+            EigenEquationSolverSettings::LinearEquationMethod EigenEquationSolverSettings::getLinearEquationSystemMethod() const {
+                std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString();
+                if (linearEquationSystemTechniqueAsString == "sparselu") {
+                    return EigenEquationSolverSettings::LinearEquationMethod::SparseLU;
+                } else if (linearEquationSystemTechniqueAsString == "bicgstab") {
+                    return EigenEquationSolverSettings::LinearEquationMethod::Bicgstab;
+                }
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected.");
+            }
+            
+            bool EigenEquationSolverSettings::isPreconditioningMethodSet() const {
+                return this->getOption(preconditionOptionName).getHasOptionBeenSet();
+            }
+            
+            EigenEquationSolverSettings::PreconditioningMethod EigenEquationSolverSettings::getPreconditioningMethod() const {
+                std::string PreconditioningMethodAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString();
+                if (PreconditioningMethodAsString == "ilu") {
+                    return EigenEquationSolverSettings::PreconditioningMethod::Ilu;
+                } else if (PreconditioningMethodAsString == "diagonal") {
+                    return EigenEquationSolverSettings::PreconditioningMethod::Diagonal;
+                } else if (PreconditioningMethodAsString == "none") {
+                    return EigenEquationSolverSettings::PreconditioningMethod::None;
+                }
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << PreconditioningMethodAsString << "' selected.");
+            }
+            
+            bool EigenEquationSolverSettings::isMaximalIterationCountSet() const {
+                return this->getOption(maximalIterationsOptionName).getHasOptionBeenSet();
+            }
+            
+            uint_fast64_t EigenEquationSolverSettings::getMaximalIterationCount() const {
+                return this->getOption(maximalIterationsOptionName).getArgumentByName("count").getValueAsUnsignedInteger();
+            }
+            
+            bool EigenEquationSolverSettings::isPrecisionSet() const {
+                return this->getOption(precisionOptionName).getHasOptionBeenSet();
+            }
+            
+            double EigenEquationSolverSettings::getPrecision() const {
+                return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble();
+            }
+            
+            bool EigenEquationSolverSettings::check() const {
+                // This list does not include the precision, because this option is shared with other modules.
+                bool optionsSet = isLinearEquationSystemMethodSet() || isPreconditioningMethodSet() || isMaximalIterationCountSet();
+                
+                STORM_LOG_WARN_COND(storm::settings::getModule<storm::settings::modules::MarkovChainSettings>().getEquationSolver() == storm::solver::EquationSolverType::Gmmxx || !optionsSet, "eigen is not selected as the preferred equation solver, so setting options for eigen might have no effect.");
+                
+                return true;
+            }
+            
+        } // namespace modules
+    } // namespace settings
+} // namespace storm

src/settings/modules/EigenEquationSolverSettings.h

@@ -0,0 +1,100 @@
+#pragma once
+
+#include "src/settings/modules/ModuleSettings.h"
+
+namespace storm {
+    namespace settings {
+        namespace modules {
+            
+            /*!
+             * This class represents the settings for Eigen.
+             */
+            class EigenEquationSolverSettings : public ModuleSettings {
+            public:
+                // An enumeration of all available methods for solving linear equations.
+                enum class LinearEquationMethod { SparseLU, Bicgstab };
+                
+                // An enumeration of all available preconditioning methods.
+                enum class PreconditioningMethod { Ilu, Diagonal, None };
+                
+                // An enumeration of all available convergence criteria.
+                enum class ConvergenceCriterion { Absolute, Relative };
+                
+                /*!
+                 * Creates a new set of Eigen settings.
+                 */
+                EigenEquationSolverSettings();
+                
+                /*!
+                 * Retrieves whether the linear equation system method has been set.
+                 *
+                 * @return True iff the linear equation system method has been set.
+                 */
+                bool isLinearEquationSystemMethodSet() const;
+                
+                /*!
+                 * Retrieves the method that is to be used for solving systems of linear equations.
+                 *
+                 * @return The method to use.
+                 */
+                LinearEquationMethod getLinearEquationSystemMethod() const;
+                
+                /*!
+                 * Retrieves whether the preconditioning method has been set.
+                 *
+                 * @return True iff the preconditioning method has been set.
+                 */
+                bool isPreconditioningMethodSet() const;
+                
+                /*!
+                 * Retrieves the method that is to be used for preconditioning solving systems of linear equations.
+                 *
+                 * @return The method to use.
+                 */
+                PreconditioningMethod getPreconditioningMethod() 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;
+                                
+                bool check() const override;
+                
+                // The name of the module.
+                static const std::string moduleName;
+                
+            private:
+                // Define the string names of the options as constants.
+                static const std::string techniqueOptionName;
+                static const std::string preconditionOptionName;
+                static const std::string maximalIterationsOptionName;
+                static const std::string maximalIterationsOptionShortName;
+                static const std::string precisionOptionName;
+            };
+            
+        } // namespace modules
+    } // namespace settings
+} // namespace storm

src/settings/modules/MarkovChainSettings.cpp

@@ -41,9 +41,9 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, engineOptionName, false, "Sets which engine is used for model building and model checking.").setShortName(engineOptionShortName)
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the engine to use. Available are {sparse, hybrid, dd, expl, abs}.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(engines)).setDefaultValueString("sparse").build()).build());
                 
-                std::vector<std::string> linearEquationSolver = {"gmm++", "native"};
+                std::vector<std::string> linearEquationSolver = {"gmm++", "native", "eigen"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, eqSolverOptionName, false, "Sets which solver is preferred for solving systems of linear equations.")
-                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the solver to prefer. Available are: gmm++ and native.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(linearEquationSolver)).setDefaultValueString("gmm++").build()).build());
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the solver to prefer. Available are: gmm++, native and eigen.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(linearEquationSolver)).setDefaultValueString("gmm++").build()).build());
                 
                 std::vector<std::string> ddLibraries = {"cudd", "sylvan"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, ddLibraryOptionName, false, "Sets which library is preferred for decision-diagram operations.")
@@ -85,6 +85,8 @@ namespace storm {
                     return storm::solver::EquationSolverType::Gmmxx;
                 } else if (equationSolverName == "native") {
                     return storm::solver::EquationSolverType::Native;
+                } else if (equationSolverName == "eigen") {
+                    return storm::solver::EquationSolverType::Eigen;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown equation solver '" << equationSolverName << "'.");
             }

src/solver/EigenLinearEquationSolver.cpp

@@ -0,0 +1,117 @@
+#include "src/solver/EigenLinearEquationSolver.h"
+
+#include "src/adapters/EigenAdapter.h"
+
+#include "src/settings/SettingsManager.h"
+#include "src/settings/modules/EigenEquationSolverSettings.h"
+
+namespace storm {
+    namespace solver {
+     
+        template<typename ValueType>
+        EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A)), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ValueType>
+        EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A)) {
+            // Get the settings object to customize linear solving.
+            storm::settings::modules::EigenEquationSolverSettings const& settings = storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>();
+            
+            // Get appropriate settings.
+            maximalNumberOfIterations = settings.getMaximalIterationCount();
+            precision = settings.getPrecision();
+            
+            // Determine the method to be used.
+            storm::settings::modules::EigenEquationSolverSettings::LinearEquationMethod methodAsSetting = settings.getLinearEquationSystemMethod();
+            if (methodAsSetting == storm::settings::modules::EigenEquationSolverSettings::LinearEquationMethod::Bicgstab) {
+                method = SolutionMethod::Bicgstab;
+            } else if (methodAsSetting == storm::settings::modules::EigenEquationSolverSettings::LinearEquationMethod::SparseLU) {
+                method = SolutionMethod::SparseLU;
+            }
+            
+            // Check which preconditioner to use.
+            storm::settings::modules::EigenEquationSolverSettings::PreconditioningMethod preconditionAsSetting = settings.getPreconditioningMethod();
+            if (preconditionAsSetting == storm::settings::modules::EigenEquationSolverSettings::PreconditioningMethod::Ilu) {
+                preconditioner = Preconditioner::Ilu;
+            } else if (preconditionAsSetting == storm::settings::modules::EigenEquationSolverSettings::PreconditioningMethod::Diagonal) {
+                preconditioner = Preconditioner::Diagonal;
+            } else if (preconditionAsSetting == storm::settings::modules::EigenEquationSolverSettings::PreconditioningMethod::None) {
+                preconditioner = Preconditioner::None;
+            }
+        }
+        
+        template<typename ValueType>
+        void EigenLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+            // Translate the vectors x and b into Eigen's format.
+            Eigen::VectorXd eigenB(b.size());
+            for (uint64_t index = 0; index < b.size(); ++index) {
+                eigenB[index] = b[index];
+            }
+            
+            Eigen::VectorXd eigenX(x.size());
+            for (uint64_t index = 0; index < x.size(); ++index) {
+                eigenX[index] = x[index];
+            }
+            
+            if (method == SolutionMethod::SparseLU) {
+                Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
+                solver.compute(*eigenA);
+                solver._solve(eigenB, eigenX);
+            } else {
+                if (preconditioner == Preconditioner::Ilu) {
+                    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::IncompleteLUT<double>> solver;
+                    solver.compute(*eigenA);
+                    solver._solve(eigenB, eigenX);
+                } else if (preconditioner == Preconditioner::Diagonal) {
+                    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::DiagonalPreconditioner<double>> solver;
+                    solver.compute(*eigenA);
+                    solver._solve(eigenB, eigenX);
+                } else {
+                    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::IdentityPreconditioner> solver;
+                    solver.compute(*eigenA);
+                    solver._solve(eigenB, eigenX);
+                }
+            }
+
+            // Translate the solution from Eigen's format into our representation.
+            for (uint64_t index = 0; index < eigenX.size(); ++index) {
+                x[index] = eigenX[index];
+            }
+        }
+        
+        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 {
+            
+            // Translate the vectors x and b into Eigen's format.
+            std::unique_ptr<Eigen::VectorXd> eigenB;
+            if (b != nullptr) {
+                eigenB = std::make_unique<Eigen::VectorXd>(b->size());
+                for (uint64_t index = 0; index < b->size(); ++index) {
+                    (*eigenB)[index] = (*b)[index];
+                }
+            }
+            Eigen::VectorXd eigenX(x.size());
+            for (uint64_t index = 0; index < x.size(); ++index) {
+                eigenX[index] = x[index];
+            }
+            
+            // Perform n matrix-vector multiplications.
+            for (uint64_t iteration = 0; iteration < n; ++iteration) {
+                eigenX = *eigenA * eigenX;
+                if (eigenB != nullptr) {
+                    eigenX += *eigenB;
+                }
+            }
+            
+            // Translate the solution from Eigen's format into our representation.
+            for (uint64_t index = 0; index < eigenX.size(); ++index) {
+                x[index] = eigenX[index];
+            }
+        }
+        
+        template class EigenLinearEquationSolver<double>;
+        
+        
+    }
+}

src/solver/EigenLinearEquationSolver.h

@@ -0,0 +1,53 @@
+#pragma once
+
+#include "src/solver/LinearEquationSolver.h"
+
+#include "src/utility/eigen.h"
+
+namespace storm {
+    namespace solver {
+        
+        /*!
+         * A class that uses the Eigen library to implement the LinearEquationSolver interface.
+         */
+        template<typename ValueType>
+        class EigenLinearEquationSolver : public LinearEquationSolver<ValueType> {
+        public:
+            enum class SolutionMethod {
+                SparseLU, Bicgstab
+            };
+            
+            enum class Preconditioner {
+                Ilu, Diagonal, None
+            };
+            
+            EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner);
+
+            EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
+            
+            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;
+            
+        private:
+            // A pointer to the original sparse matrix given to this solver.
+            storm::storage::SparseMatrix<ValueType> const* originalA;
+            
+            // The (eigen) matrix associated with this equation solver.
+            std::unique_ptr<Eigen::SparseMatrix<double>> eigenA;
+
+            // The method to use for solving linear equation systems.
+            SolutionMethod method;
+            
+            // The required precision for the iterative methods.
+            double precision;
+            
+            // 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;
+        };
+        
+    }
+}

src/solver/GmmxxLinearEquationSolver.h

@@ -97,7 +97,7 @@ namespace storm {
             // The preconditioner to use when solving the linear equation system.
             Preconditioner preconditioner;
 
-            // Sets whether the relative or absolute error is to be considered for convergence detection. Not that this
+            // Sets whether the relative or absolute error is to be considered for convergence detection. Note that this
             // only applies to the Jacobi method for this solver.
             bool relative;
 

src/solver/SolverSelectionOptions.h

@@ -9,7 +9,7 @@ namespace storm {
         ExtendEnumsWithSelectionField(MinMaxTechnique, PolicyIteration, ValueIteration)
         
         ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk)
-        ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Topological)
+        ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Topological)
         ExtendEnumsWithSelectionField(SmtSolverType, Z3, Mathsat)
     }
 } 

src/utility/ConversionHelper.cpp

@@ -1,5 +0,0 @@
-#include "ConversionHelper.h"
-
-std::vector<unsigned long long, std::allocator<unsigned long long>>* storm::utility::ConversionHelper::toUnsignedLongLong(std::vector<uint_fast64_t, std::allocator<uint_fast64_t>>* vectorPtr) {
-	return reinterpret_cast<std::vector<unsigned long long, std::allocator<unsigned long long>> *>(vectorPtr);
-}

src/utility/ConversionHelper.h

@@ -1,38 +0,0 @@
-/*
- * ConversionHelper.h
- *
- *  Created on: 14.09.2013
- *      Author: Philipp Berger
- * 
- * WARNING: This file REQUIRES -no-strict-aliasing!
- */
-
-#ifndef STORM_UTILITY_CONVERSIONHELPER_H_
-#define STORM_UTILITY_CONVERSIONHELPER_H_
-
-#include <iostream>
-#include <vector>
-#include <cstdint>
-
-static_assert(sizeof(unsigned long long) == sizeof(uint_fast64_t), "This program uses the GMM Backend and therefor requires unsigned long long and uint_fast64_t to be of the same size!");
-
-
-namespace storm {
-	namespace utility {
-
-		class ConversionHelper {
-		public:
-			/*!
-			* Converts a pointer to a std::vector<uint_fast64_t> to std::vector<unsigned long long>
-			*/
-			static std::vector<unsigned long long, std::allocator<unsigned long long>>* toUnsignedLongLong(std::vector<uint_fast64_t, std::allocator<uint_fast64_t>>* vectorPtr);
-
-		private:
-			ConversionHelper() {}
-			ConversionHelper(ConversionHelper& other) {}
-			~ConversionHelper() {}
-		};
-	}
-}
-
-#endif // STORM_UTILITY_CONVERSIONHELPER_H_

src/utility/eigen.h

@@ -0,0 +1,3 @@
+#pragma once
+
+#include <Eigen/Sparse>

src/utility/solver.cpp

@@ -7,6 +7,7 @@
 #include "src/solver/SymbolicGameSolver.h"
 #include "src/solver/NativeLinearEquationSolver.h"
 #include "src/solver/GmmxxLinearEquationSolver.h"
+#include "src/solver/EigenLinearEquationSolver.h"
 #include "src/solver/GameSolver.h"
 
 #include "src/solver/NativeMinMaxLinearEquationSolver.h"
@@ -50,6 +51,7 @@ namespace storm {
                 switch (equationSolver) {
                     case storm::solver::EquationSolverType::Gmmxx: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
                     case storm::solver::EquationSolverType::Native: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix));
+                    case storm::solver::EquationSolverType::Eigen: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::EigenLinearEquationSolver<ValueType>(matrix));
                     default: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
                 }
             }
@@ -59,6 +61,11 @@ namespace storm {
                 return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
             }
 
+            template<typename ValueType>
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EigenLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+                return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::EigenLinearEquationSolver<ValueType>(matrix));
+            }
+            
             template<typename ValueType>
             NativeLinearEquationSolverFactory<ValueType>::NativeLinearEquationSolverFactory() {
                 switch (storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getLinearEquationSystemMethod()) {

src/utility/solver.h

@@ -106,6 +106,12 @@ namespace storm {
                 virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
             };
             
+            template<typename ValueType>
+            class EigenLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
+            public:
+                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            };
+            
             template<typename ValueType>
             class MinMaxLinearEquationSolverFactory {
             public:

test/functional/solver/EigenLinearEquationSolverTest.cpp

@@ -0,0 +1,170 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/solver/EigenLinearEquationSolver.h"
+#include "src/settings/SettingsManager.h"
+
+#include "src/settings/modules/GmmxxEquationSolverSettings.h"
+
+TEST(EigenLinearEquationSolver, SolveWithStandardOptions) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, 2));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 4));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, -2));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, 4));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, 3));
+    
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(3);
+    std::vector<double> b = {16, -4, -7};
+    
+    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<double> solver(A));
+    
+    storm::solver::EigenLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+}
+
+TEST(EigenLinearEquationSolver, SparseLU) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, 2));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 4));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, -2));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, 4));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, 3));
+    
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(3);
+    std::vector<double> b = {16, -4, -7};
+    
+    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::SparseLU, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None));
+    
+    storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::SparseLU, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+}
+
+TEST(EigenLinearEquationSolver, BiCGSTAB) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, 2));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 4));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, -2));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, 4));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, 3));
+    
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(3);
+    std::vector<double> b = {16, -4, -7};
+    
+    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None));
+    
+    storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+}
+
+TEST(EigenLinearEquationSolver, BiCGSTAB_Ilu) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, 2));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 4));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, -2));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, 4));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, 3));
+    
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(3);
+    std::vector<double> b = {16, -4, -7};
+    
+    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Ilu));
+    
+    storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Ilu);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+}
+
+TEST(EigenLinearEquationSolver, BiCGSTAB_Diagonal) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 0, 2));
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 4));
+    ASSERT_NO_THROW(builder.addNextValue(0, 2, -2));
+    ASSERT_NO_THROW(builder.addNextValue(1, 0, 4));
+    ASSERT_NO_THROW(builder.addNextValue(1, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 0, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 1, -1));
+    ASSERT_NO_THROW(builder.addNextValue(2, 2, 3));
+    
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(3);
+    std::vector<double> b = {16, -4, -7};
+    
+    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Diagonal));
+    
+    storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Diagonal);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+    ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+}
+
+TEST(EigenLinearEquationSolver, MatrixVectorMultplication) {
+    ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
+    storm::storage::SparseMatrixBuilder<double> builder;
+    ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(0, 4, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(1, 2, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(1, 4, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 3, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(2, 4, 0.5));
+    ASSERT_NO_THROW(builder.addNextValue(3, 4, 1));
+    ASSERT_NO_THROW(builder.addNextValue(4, 4, 1));
+    
+    storm::storage::SparseMatrix<double> A;
+    ASSERT_NO_THROW(A = builder.build());
+    
+    std::vector<double> x(5);
+    x[4] = 1;
+    
+    storm::solver::EigenLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(x, nullptr, 4));
+    ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
+}

test/functional/solver/GmmxxLinearEquationSolverTest.cpp

@@ -167,7 +167,7 @@ TEST(GmmxxLinearEquationSolver, gmresilu) {
     
     ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Ilu, true, 50));
     
-    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50);
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Ilu, true, 50);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
@@ -195,7 +195,7 @@ TEST(GmmxxLinearEquationSolver, gmresdiag) {
     
     ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Diagonal, true, 50));
     
-    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50);
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Diagonal, true, 50);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());