From b67e3d6e7bd2d17ca296d0480b0b53c07a7d6b56 Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Fri, 5 Aug 2016 11:27:36 +0200
Subject: [PATCH] added 'convergence' (rather success) checks for Eigen solver

Former-commit-id: 25a6fb1d77945935fe43153ee9e1e0e53261bc5e
---
 src/solver/EigenLinearEquationSolver.cpp | 18 +++++++++++++++---
 src/solver/EigenLinearEquationSolver.h   |  2 +-
 2 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp
index 47117b9a3..eb19b57e2 100644
--- a/src/solver/EigenLinearEquationSolver.cpp
+++ b/src/solver/EigenLinearEquationSolver.cpp
@@ -127,7 +127,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void EigenLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+        bool EigenLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             // Map the input vectors to Eigen's format.
             auto eigenX = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(x.data(), x.size());
             auto eigenB = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b.data(), b.size());
@@ -146,18 +146,21 @@ namespace storm {
                         solver.setTolerance(this->getSettings().getPrecision());
                         solver.setMaxIterations(this->getSettings().getMaximalNumberOfIterations());
                         eigenX = solver.solveWithGuess(eigenB, eigenX);
+                        return solver.info() == Eigen::ComputationInfo::Success;
                     } else if (preconditioner == EigenLinearEquationSolverSettings<ValueType>::Preconditioner::Diagonal) {
                         Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::DiagonalPreconditioner<ValueType>> solver;
                         solver.setTolerance(this->getSettings().getPrecision());
                         solver.setMaxIterations(this->getSettings().getMaximalNumberOfIterations());
                         solver.compute(*this->eigenA);
                         eigenX = solver.solveWithGuess(eigenB, eigenX);
+                        return solver.info() == Eigen::ComputationInfo::Success;
                     } else {
                         Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::IdentityPreconditioner> solver;
                         solver.setTolerance(this->getSettings().getPrecision());
                         solver.setMaxIterations(this->getSettings().getMaximalNumberOfIterations());
                         solver.compute(*this->eigenA);
                         eigenX = solver.solveWithGuess(eigenB, eigenX);
+                        return solver.info() == Eigen::ComputationInfo::Success;
                     }
                 } else if (solutionMethod == EigenLinearEquationSolverSettings<ValueType>::SolutionMethod::DGMRES) {
                     if (preconditioner == EigenLinearEquationSolverSettings<ValueType>::Preconditioner::Ilu) {
@@ -167,6 +170,7 @@ namespace storm {
                         solver.set_restart(this->getSettings().getNumberOfIterationsUntilRestart());
                         solver.compute(*this->eigenA);
                         eigenX = solver.solveWithGuess(eigenB, eigenX);
+                        return solver.info() == Eigen::ComputationInfo::Success;
                     } else if (preconditioner == EigenLinearEquationSolverSettings<ValueType>::Preconditioner::Diagonal) {
                         Eigen::DGMRES<Eigen::SparseMatrix<ValueType>, Eigen::DiagonalPreconditioner<ValueType>> solver;
                         solver.setTolerance(this->getSettings().getPrecision());
@@ -174,6 +178,7 @@ namespace storm {
                         solver.set_restart(this->getSettings().getNumberOfIterationsUntilRestart());
                         solver.compute(*this->eigenA);
                         eigenX = solver.solveWithGuess(eigenB, eigenX);
+                        return solver.info() == Eigen::ComputationInfo::Success;
                     } else {
                         Eigen::DGMRES<Eigen::SparseMatrix<ValueType>, Eigen::IdentityPreconditioner> solver;
                         solver.setTolerance(this->getSettings().getPrecision());
@@ -181,6 +186,7 @@ namespace storm {
                         solver.set_restart(this->getSettings().getNumberOfIterationsUntilRestart());
                         solver.compute(*this->eigenA);
                         eigenX = solver.solveWithGuess(eigenB, eigenX);
+                        return solver.info() == Eigen::ComputationInfo::Success;
                     }
                 } else if (solutionMethod == EigenLinearEquationSolverSettings<ValueType>::SolutionMethod::GMRES) {
                     if (preconditioner == EigenLinearEquationSolverSettings<ValueType>::Preconditioner::Ilu) {
@@ -190,6 +196,7 @@ namespace storm {
                         solver.set_restart(this->getSettings().getNumberOfIterationsUntilRestart());
                         solver.compute(*this->eigenA);
                         eigenX = solver.solveWithGuess(eigenB, eigenX);
+                        return solver.info() == Eigen::ComputationInfo::Success;
                     } else if (preconditioner == EigenLinearEquationSolverSettings<ValueType>::Preconditioner::Diagonal) {
                         Eigen::GMRES<Eigen::SparseMatrix<ValueType>, Eigen::DiagonalPreconditioner<ValueType>> solver;
                         solver.setTolerance(this->getSettings().getPrecision());
@@ -197,6 +204,7 @@ namespace storm {
                         solver.set_restart(this->getSettings().getNumberOfIterationsUntilRestart());
                         solver.compute(*this->eigenA);
                         eigenX = solver.solveWithGuess(eigenB, eigenX);
+                        return solver.info() == Eigen::ComputationInfo::Success;
                     } else {
                         Eigen::GMRES<Eigen::SparseMatrix<ValueType>, Eigen::IdentityPreconditioner> solver;
                         solver.setTolerance(this->getSettings().getPrecision());
@@ -204,9 +212,11 @@ namespace storm {
                         solver.set_restart(this->getSettings().getNumberOfIterationsUntilRestart());
                         solver.compute(*this->eigenA);
                         eigenX = solver.solveWithGuess(eigenB, eigenX);
+                        return solver.info() == Eigen::ComputationInfo::Success;
                     }
                 }
             }
+            return false;
         }
         
         template<typename ValueType>
@@ -260,7 +270,7 @@ namespace storm {
 #ifdef STORM_HAVE_CARL
         // Specialization for storm::RationalNumber
         template<>
-        void EigenLinearEquationSolver<storm::RationalNumber>::solveEquations(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b) const {
+        bool EigenLinearEquationSolver<storm::RationalNumber>::solveEquations(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b) const {
             // Map the input vectors to Eigen's format.
             auto eigenX = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(x.data(), x.size());
             auto eigenB = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b.data(), b.size());
@@ -268,11 +278,12 @@ namespace storm {
             Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalNumber>, Eigen::COLAMDOrdering<int>> solver;
             solver.compute(*eigenA);
             solver._solve_impl(eigenB, eigenX);
+            return solver.info() == Eigen::ComputationInfo::Success;
         }
         
         // Specialization for storm::RationalFunction
         template<>
-        void EigenLinearEquationSolver<storm::RationalFunction>::solveEquations(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b) const {
+        bool EigenLinearEquationSolver<storm::RationalFunction>::solveEquations(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b) const {
             // Map the input vectors to Eigen's format.
             auto eigenX = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(x.data(), x.size());
             auto eigenB = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b.data(), b.size());
@@ -280,6 +291,7 @@ namespace storm {
             Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalFunction>, Eigen::COLAMDOrdering<int>> solver;
             solver.compute(*eigenA);
             solver._solve_impl(eigenB, eigenX);
+            return solver.info() == Eigen::ComputationInfo::Success;
         }
 #endif
 
diff --git a/src/solver/EigenLinearEquationSolver.h b/src/solver/EigenLinearEquationSolver.h
index 688dcae64..ee20f5a30 100644
--- a/src/solver/EigenLinearEquationSolver.h
+++ b/src/solver/EigenLinearEquationSolver.h
@@ -66,7 +66,7 @@ namespace storm {
             virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
             virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
             
-            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
+            virtual bool solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
             virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
 
             EigenLinearEquationSolverSettings<ValueType>& getSettings();