Extended interface of linEqSolvers a little,

eq solver is now called repeatedly with increased precision, until the result is good enough..


Former-commit-id: b076950dc8

src/solver/GmmxxLinearEquationSolver.cpp

@@ -55,7 +55,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void GmmxxLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+        bool GmmxxLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner '" << preconditionerToString() << "' (max. " << maximalNumberOfIterations << " iterations).");
             if (method == SolutionMethod::Jacobi && preconditioner != Preconditioner::None) {
                 LOG4CPLUS_WARN(logger, "Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored.");
@@ -94,8 +94,10 @@ namespace storm {
                 // Check if the solver converged and issue a warning otherwise.
                 if (iter.converged()) {
                     LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
+                    return true;
                 } else {
-                    LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
+                    LOG4CPLUS_INFO(logger, "Iterative solver did not converge.");
+                    return false;
                 }
             } else if (method == SolutionMethod::Jacobi) {
                 uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(*originalA, x, b, multiplyResult);
@@ -103,9 +105,14 @@ namespace storm {
                 // Check if the solver converged and issue a warning otherwise.
                 if (iterations < maximalNumberOfIterations) {
                     LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
+                    return true;
                 } else {
-                    LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
+                    LOG4CPLUS_INFO(logger, "Iterative solver did not converge.");
+                    return false;
                 }
+            } else {
+                STORM_LOG_ERROR("The selected method is not supported.");
+                return false;
             }
         }
         

src/solver/GmmxxLinearEquationSolver.h

@@ -46,10 +46,22 @@ namespace storm {
              */
             GmmxxLinearEquationSolver(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 bool solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
+            virtual void setPrecision(double precision) override{
+                this->precision = precision;
+            }
+            
+            virtual void setIterations(uint_fast64_t maximalNumberOfIterations) override{
+                this->maximalNumberOfIterations = maximalNumberOfIterations;
+            }
+            
+            virtual void setRelative(bool relative) override{
+                this->relative = relative;
+            }
+            
         private:
             /*!
              * Solves the linear equation system A*x = b given by the parameters using the Jacobi method.

src/solver/LinearEquationSolver.h

@@ -29,8 +29,10 @@ namespace storm {
              * @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.
+             * 
+             * @return true iff solving was successful (e.g. iterative methods converged)
              */
-            virtual void solveEquationSystem(std::vector<Type>& x, std::vector<Type> const& b, std::vector<Type>* multiplyResult = nullptr) const = 0;
+            virtual bool solveEquationSystem(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
@@ -45,6 +47,12 @@ namespace storm {
              * vector must be equal to the number of rows of A.
              */
             virtual void performMatrixVectorMultiplication(std::vector<Type>& x, std::vector<Type> const* b = nullptr, uint_fast64_t n = 1, std::vector<Type>* multiplyResult = nullptr) const = 0;
+            
+            virtual void setPrecision(double precision) = 0;
+            virtual void setIterations(uint_fast64_t maximalNumberOfIterations) = 0;
+            
+            virtual void setRelative(bool relative) = 0;
+            
         };
         
     } // namespace solver

src/solver/NativeLinearEquationSolver.cpp

@@ -28,7 +28,7 @@ namespace storm {
         }
                 
         template<typename ValueType>
-        void NativeLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+        bool NativeLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             if (method == NativeLinearEquationSolverSolutionMethod::SOR || method == NativeLinearEquationSolverSolutionMethod::GaussSeidel) {
                 // Define the omega used for SOR.
                 ValueType omega = method == NativeLinearEquationSolverSolutionMethod::SOR ? storm::settings::nativeEquationSolverSettings().getOmega() : storm::utility::one<ValueType>();
@@ -66,6 +66,9 @@ namespace storm {
                 if (!tmpXProvided) {
                     delete tmpX;
                 }
+                
+                return converged;
+                
             } else {
                 // Get a Jacobi decomposition of the matrix A.
                 std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
@@ -113,6 +116,8 @@ namespace storm {
                 if (!multiplyResultProvided) {
                     delete copyX;
                 }
+                
+                return converged;
             }
         }
         

src/solver/NativeLinearEquationSolver.h

@@ -38,10 +38,23 @@ namespace storm {
              */
             NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
                         
-            virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            virtual bool solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
+            virtual void setPrecision(double precision) override{
+                this->precision = precision;
+            }
+            
+            virtual void setIterations(uint_fast64_t maximalNumberOfIterations) override{
+                this->maximalNumberOfIterations = maximalNumberOfIterations;
+            }
+            
+            virtual void setRelative(bool relative) override{
+                this->relative = relative;
+            }
+            
+
         private:
             /*!
              * Retrieves the string representation of the solution method associated with this solver.

src/storage/SparseMatrix.cpp

@@ -268,7 +268,6 @@ namespace storm {
         template<typename ValueType>
         SparseMatrix<ValueType>::SparseMatrix(SparseMatrix<ValueType> const& other) : rowCount(other.rowCount), columnCount(other.columnCount), entryCount(other.entryCount), nonzeroEntryCount(other.nonzeroEntryCount), columnsAndValues(other.columnsAndValues), rowIndications(other.rowIndications), nontrivialRowGrouping(other.nontrivialRowGrouping), rowGroupIndices(other.rowGroupIndices) {
             // Intentionally left empty.
-            std::cout << "!!!! copying matrix (constructor)" << std::endl;
         }
         
         template<typename ValueType>
@@ -276,7 +275,6 @@ namespace storm {
             storm::storage::BitVector rowConstraint(other.getRowCount(), true);
             storm::storage::BitVector columnConstraint(other.getColumnCount(), true);
             *this = other.getSubmatrix(false, rowConstraint, columnConstraint, insertDiagonalElements);
-            std::cout << "!!!! copying matrix (insertingDiagElements)" << std::endl;
         }
         
         template<typename ValueType>
@@ -290,7 +288,6 @@ namespace storm {
         template<typename ValueType>
         SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, ValueType>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(rowGroupIndices) {
             this->updateNonzeroEntryCount();
-            std::cout << "!!!! copying matrix (copying Ingredients)" << std::endl;
         }
         
         template<typename ValueType>
@@ -300,7 +297,6 @@ namespace storm {
         
         template<typename ValueType>
         SparseMatrix<ValueType>& SparseMatrix<ValueType>::operator=(SparseMatrix<ValueType> const& other) {
-            std::cout << "!!!! copying matrix (operator=)" << std::endl;
             // Only perform assignment if source and target are not the same.
             if (this != &other) {
                 rowCount = other.rowCount;

src/utility/policyguessing.cpp

@@ -201,25 +201,30 @@ namespace storm {
                                         ){
                 //Get the submatrix/subvector A,x, and b and invoke linear equation solver
                 storm::storage::SparseMatrix<ValueType> subA = A.getSubmatrix(true, probGreater0States, probGreater0States, true);
-                subA.convertToEquationSystem();
+                storm::storage::SparseMatrix<ValueType> eqSysA(subA);
+                eqSysA.convertToEquationSystem();
                 std::vector<ValueType> subX(probGreater0States.getNumberOfSetBits());
                 storm::utility::vector::selectVectorValues(subX, probGreater0States, x);
                 std::vector<ValueType> subB(probGreater0States.getNumberOfSetBits());
                 storm::utility::vector::selectVectorValues(subB, probGreater0States, b);
-                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory<ValueType>().create(subA);
-                linEqSysSolver->solveEquationSystem(subX, subB);
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory<ValueType>().create(eqSysA);
+                linEqSysSolver->setRelative(relative);
                 
-                //Performing a couple of iterations makes the result "stable" when doing value iteration, i.e.,
-                //if the given equation system is induced by an optimal policy, value iteration will terminate after the first iteration.
-                subA.convertToEquationSystem();
-                linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory<ValueType>().create(subA);
+                std::size_t iterations = 0;
                 std::vector<ValueType> copyX(subX.size());
-                std::size_t iterations = 100000; //don't run into an endless loop...
+                ValueType newPrecision = precision;
                 do {
-                    linEqSysSolver->performMatrixVectorMultiplication(subX, &subB, 10, &copyX);
-                    --iterations;
-                } while(!storm::utility::vector::equalModuloPrecision(subX, copyX, precision, relative) && iterations>0);
-                STORM_LOG_WARN_COND(iterations>0, "Iterations on result of linear equation solver did not converge.");
+                    linEqSysSolver->setPrecision(newPrecision);
+                    if(!linEqSysSolver->solveEquationSystem(subX, subB)){
+                        break; //Solver did not converge.. so we have to go on with the current solution.
+                    }
+                    subA.multiplyWithVector(subX,copyX);
+                    storm::utility::vector::addVectors(copyX, subB, copyX); // = Ax + b
+                    ++iterations;
+                    newPrecision *= 0.1;
+                } while(!storm::utility::vector::equalModuloPrecision(subX, copyX, precision, relative) && iterations<30);
+                
+                STORM_LOG_DEBUG("Required to increase the precision " << iterations << " times in order to obtain a precise result");
                 //fill in the result
                 storm::utility::vector::setVectorValues(x, probGreater0States, subX);
                 storm::utility::vector::setVectorValues(x, (~probGreater0States), prob0Value);

src/utility/policyguessing.h

@@ -150,11 +150,16 @@ namespace storm {
              * However, actual target states are already filtered out.
              * To ensure a unique solution, we also need to filter out the "prob0"-states.
              * 
+             * If the result does not satisfy "Ax+b = x (modulo precision)", the solver is executed
+             * again with increased precision.
+             *
              * @param A the matrix of the equation system
              * @param x The initial guess of the solution.
              * @param b The vector of the equation system
              * @param targetChoices marks the rows in the matrix that have a positive probability to lead to a target state
              * @param prob0Value the value that is assigned to the states that have probability zero to reach a target
+             * @param const& precision The precision to be used by the solver
+             * @param relative sets whether to consider relative errors
              * @return The solution vector in the form of the vector x.
              */
             template<typename ValueType>