First steps toward exact linear equation solver

Former-commit-id: 669af13b84

src/adapters/CarlAdapter.h

@@ -4,6 +4,8 @@
 // Include config to know whether CARL is available or not.
 #include "storm-config.h"
 
+#include <boost/multiprecision/gmp.hpp>
+
 #ifdef STORM_HAVE_CARL
 
 #include <carl/numbers/numbers.h>
@@ -42,9 +44,12 @@ namespace carl {
 }
 
 namespace storm {
-    typedef cln::cl_RA RationalNumber;
+//    typedef boost::multiprecision::gmp_rational RationalNumber;
+    typedef mpq_class RationalNumber;
+    
+    typedef cln::cl_RA CarlRationalNumber;
     typedef carl::Variable Variable;
-    typedef carl::MultivariatePolynomial<RationalNumber> RawPolynomial;
+    typedef carl::MultivariatePolynomial<CarlRationalNumber> RawPolynomial;
     typedef carl::FactorizedPolynomial<RawPolynomial> Polynomial;
     typedef carl::Relation CompareRelation;
     

src/solver/EliminationLinearEquationSolver.cpp

@@ -0,0 +1,62 @@
+#include "src/solver/EliminationLinearEquationSolver.h"
+
+#include "src/utility/vector.h"
+
+namespace storm {
+    namespace solver {
+        template<typename ValueType>
+        EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : A(A) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ValueType>
+        void EliminationLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+            // TODO: implement state-elimination here.
+        }
+        
+        template<typename ValueType>
+        void EliminationLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+            // Set up some temporary variables so that we can just swap pointers instead of copying the result after
+            // each iteration.
+            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) {
+                A.multiplyWithVector(*currentX, *nextX);
+                std::swap(nextX, currentX);
+                
+                // If requested, add an offset to the current result vector.
+                if (b != nullptr) {
+                    storm::utility::vector::addVectors(*currentX, *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 class EliminationLinearEquationSolver<double>;
+        
+        // TODO: make this work with the proper implementation of solveEquationSystem.
+        template class EliminationLinearEquationSolver<storm::RationalNumber>;
+        template class EliminationLinearEquationSolver<storm::RationalFunction>;
+        
+    }
+}
+

src/solver/EliminationLinearEquationSolver.h

@@ -0,0 +1,34 @@
+#ifndef STORM_SOLVER_ELIMINATIONLINEAREQUATIONSOLVER_H_
+#define STORM_SOLVER_ELIMINATIONLINEAREQUATIONSOLVER_H_
+
+#include "src/solver/LinearEquationSolver.h"
+
+namespace storm {
+    namespace solver {
+        /*!
+         * A class that uses gaussian elimination to implement the LinearEquationSolver interface. In particular
+         */
+        template<typename ValueType>
+        class EliminationLinearEquationSolver : public LinearEquationSolver<ValueType> {
+        public:
+            
+            /*!
+             * Constructs a linear equation solver.
+             *
+             * @param A The matrix defining the coefficients of the linear equation system.
+             */
+            EliminationLinearEquationSolver(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 reference to the original matrix used for this equation solver.
+            storm::storage::SparseMatrix<ValueType> const& A;
+            
+        };
+    }
+}
+
+#endif /* STORM_SOLVER_ELIMINATIONLINEAREQUATIONSOLVER_H_ */

src/solver/NativeLinearEquationSolver.h

@@ -10,6 +10,7 @@ namespace storm {
         enum class NativeLinearEquationSolverSolutionMethod {
             Jacobi, GaussSeidel, SOR
         };
+        
         /*!
          * A class that uses StoRM's native matrix operations to implement the LinearEquationSolver interface.
          */
@@ -17,7 +18,6 @@ namespace storm {
         class NativeLinearEquationSolver : public LinearEquationSolver<ValueType> {
         public:
             
-            
             /*!
              * Constructs a linear equation solver with parameters being set according to the settings object.
              *

src/storage/FlexibleSparseMatrix.cpp

@@ -0,0 +1,88 @@
+#include "src/storage/FlexibleSparseMatrix.h"
+
+#include "src/storage/SparseMatrix.h"
+#include "src/storage/BitVector.h"
+
+namespace storm {
+    namespace storage {
+        template<typename ValueType>
+        FlexibleSparseMatrix<ValueType>::FlexibleSparseMatrix(index_type rows) : data(rows) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ValueType>
+        void FlexibleSparseMatrix<ValueType>::reserveInRow(index_type row, index_type numberOfElements) {
+            this->data[row].reserve(numberOfElements);
+        }
+        
+        template<typename ValueType>
+        typename FlexibleSparseMatrix<ValueType>::row_type& FlexibleSparseMatrix<ValueType>::getRow(index_type index) {
+            return this->data[index];
+        }
+        
+        template<typename ValueType>
+        typename FlexibleSparseMatrix<ValueType>::row_type const& FlexibleSparseMatrix<ValueType>::getRow(index_type index) const {
+            return this->data[index];
+        }
+        
+        template<typename ValueType>
+        typename FlexibleSparseMatrix<ValueType>::index_type FlexibleSparseMatrix<ValueType>::getNumberOfRows() const {
+            return this->data.size();
+        }
+        
+        template<typename ValueType>
+        bool FlexibleSparseMatrix<ValueType>::hasSelfLoop(storm::storage::sparse::state_type state) {
+            for (auto const& entry : this->getRow(state)) {
+                if (entry.getColumn() < state) {
+                    continue;
+                } else if (entry.getColumn() > state) {
+                    return false;
+                } else if (entry.getColumn() == state) {
+                    return true;
+                }
+            }
+            return false;
+        }
+        
+        template<typename ValueType>
+        void FlexibleSparseMatrix<ValueType>::print() const {
+            for (uint_fast64_t index = 0; index < this->data.size(); ++index) {
+                std::cout << index << " - ";
+                for (auto const& element : this->getRow(index)) {
+                    std::cout << "(" << element.getColumn() << ", " << element.getValue() << ") ";
+                }
+                std::cout << std::endl;
+            }
+        }
+        
+        template<typename ValueType>
+        bool FlexibleSparseMatrix<ValueType>::empty() const {
+            for (auto const& row : this->data) {
+                if (!row.empty()) {
+                    return false;
+                }
+            }
+            return true;
+        }
+        
+        template<typename ValueType>
+        void FlexibleSparseMatrix<ValueType>::filter(storm::storage::BitVector const& rowFilter, storm::storage::BitVector const& columnFilter) {
+            for (uint_fast64_t rowIndex = 0; rowIndex < this->data.size(); ++rowIndex) {
+                auto& row = this->data[rowIndex];
+                if (!rowFilter.get(rowIndex)) {
+                    row.clear();
+                    row.shrink_to_fit();
+                    continue;
+                }
+                row_type newRow;
+                for (auto const& element : row) {
+                    if (columnFilter.get(element.getColumn())) {
+                        newRow.push_back(element);
+                    }
+                }
+                row = std::move(newRow);
+            }
+        }
+
+    }
+}

src/storage/FlexibleSparseMatrix.h

@@ -0,0 +1,65 @@
+#ifndef STORM_STORAGE_FLEXIBLESPARSEMATRIX_H_
+#define STORM_STORAGE_FLEXIBLESPARSEMATRIX_H_
+
+#include <cstdint>
+#include <vector>
+
+#include "src/storage/sparse/StateType.h"
+
+namespace storm {
+    namespace storage {
+        template <typename IndexType, typename ValueType>
+        class MatrixEntry;
+        
+        class BitVector;
+        
+        template<typename ValueType>
+        class FlexibleSparseMatrix {
+        public:
+            // TODO: make this class a bit more consistent with the big sparse matrix and improve it:
+            // * rename getNumberOfRows -> getRowCount
+            // * store number of columns to also provide getColumnCount
+            // * rename hasSelfLoop -> rowHasDiagonalElement
+            // * add output iterator and improve the way the matrix is printed
+            // * add conversion functionality from/to sparse matrix
+            // * add documentation
+            // * rename filter to something more appropriate (getSubmatrix?)
+            // * add stuff like clearRow, multiplyRowWithScalar
+            
+            typedef uint_fast64_t index_type;
+            typedef ValueType value_type;
+            typedef std::vector<storm::storage::MatrixEntry<index_type, value_type>> row_type;
+            typedef typename row_type::iterator iterator;
+            typedef typename row_type::const_iterator const_iterator;
+            
+            FlexibleSparseMatrix() = default;
+            FlexibleSparseMatrix(index_type rows);
+            
+            void reserveInRow(index_type row, index_type numberOfElements);
+            
+            row_type& getRow(index_type);
+            row_type const& getRow(index_type) const;
+            
+            index_type getNumberOfRows() const;
+            
+            void print() const;
+            
+            bool empty() const;
+            
+            void filter(storm::storage::BitVector const& rowFilter, storm::storage::BitVector const& columnFilter);
+            
+            /*!
+             * Checks whether the given state has a self-loop with an arbitrary probability in the probability matrix.
+             *
+             * @param state The state for which to check whether it possesses a self-loop.
+             * @return True iff the given state has a self-loop with an arbitrary probability in the probability matrix.
+             */
+            bool hasSelfLoop(storm::storage::sparse::state_type state);
+            
+        private:
+            std::vector<row_type> data;
+        };
+    }
+}
+
+#endif /* STORM_STORAGE_FLEXIBLESPARSEMATRIX_H_ */

src/storage/expressions/ToRationalFunctionVisitor.cpp

@@ -88,12 +88,12 @@ namespace storm {
         
         template<typename RationalFunctionType>
         boost::any ToRationalFunctionVisitor<RationalFunctionType>::visit(IntegerLiteralExpression const& expression) {
-            return RationalFunctionType(carl::rationalize<storm::RationalNumber>(static_cast<size_t>(expression.getValue())));
+            return RationalFunctionType(carl::rationalize<storm::CarlRationalNumber>(static_cast<size_t>(expression.getValue())));
         }
         
         template<typename RationalFunctionType>
         boost::any ToRationalFunctionVisitor<RationalFunctionType>::visit(DoubleLiteralExpression const& expression) {
-            return RationalFunctionType(carl::rationalize<storm::RationalNumber>(expression.getValue()));
+            return RationalFunctionType(carl::rationalize<storm::CarlRationalNumber>(expression.getValue()));
         }
 
         template class ToRationalFunctionVisitor<storm::RationalFunction>;

src/utility/vector.h

@@ -269,7 +269,7 @@ namespace storm {
              * @param target The target vector.
              */
             template<class InValueType1, class InValueType2, class OutValueType>
-            void applyPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target, std::function<OutValueType (InValueType1 const&, InValueType2 const&)> function) {
+            void applyPointwise(std::vector<InValueType1> const& firstOperand, std::vector<InValueType2> const& secondOperand, std::vector<OutValueType>& target, std::function<OutValueType (InValueType1 const&, InValueType2 const&)> const& function) {
 #ifdef STORM_HAVE_INTELTBB
                 tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
                                   [&](tbb::blocked_range<uint_fast64_t> const& range) {