From 4e14ecb869cbf715591569c62246bc84c9afd4a8 Mon Sep 17 00:00:00 2001 From: dehnert Date: Sun, 19 Jun 2016 18:58:31 +0200 Subject: [PATCH] made elimination-based linear solver work in an alpha version. changed minor things in Eigen's SparseLU implementation to make it work with rational numbers and rational functions Former-commit-id: e5622bd98116836ffb2bfd09b7b22f1185f25960 --- .../Eigen/src/SparseCore/SparseMatrix.h | 4 +- .../eigen-3.2.6/Eigen/src/SparseLU/SparseLU.h | 2 +- .../Eigen/src/SparseLU/SparseLU_column_bmod.h | 2 +- .../src/SparseLU/SparseLU_copy_to_ucol.h | 2 +- .../Eigen/src/SparseLU/SparseLU_panel_bmod.h | 6 +- .../Eigen/src/SparseLU/SparseLU_pivotL.h | 47 +-- src/adapters/EigenAdapter.cpp | 2 + .../SparseDtmcEliminationModelChecker.cpp | 6 +- src/settings/modules/MarkovChainSettings.cpp | 6 +- src/solver/EigenLinearEquationSolver.cpp | 138 +++++++- src/solver/EigenLinearEquationSolver.h | 49 ++- .../EliminationLinearEquationSolver.cpp | 61 +++- src/solver/SolverSelectionOptions.h | 2 +- ...tor.cpp => ConditionalStateEliminator.cpp} | 31 +- ...minator.h => ConditionalStateEliminator.h} | 11 +- .../PrioritizedStateEliminator.cpp | 3 +- .../stateelimination/StateEliminator.cpp | 1 + src/storage/FlexibleSparseMatrix.cpp | 36 +- src/storage/FlexibleSparseMatrix.h | 4 +- src/utility/constants.cpp | 28 +- src/utility/solver.cpp | 15 +- src/utility/solver.h | 16 +- src/utility/stateelimination.cpp | 7 +- src/utility/stateelimination.h | 3 +- .../EigenDtmcPrctlModelCheckerTest.cpp | 321 ++++++++++++++++++ .../EigenDtmcPrctlModelCheckerTest.hpp | 14 + .../solver/EigenLinearEquationSolverTest.cpp | 91 ++++- .../EliminationLinearEquationSolverTest.cpp | 58 ++++ 28 files changed, 849 insertions(+), 117 deletions(-) rename src/solver/stateelimination/{ConditionalEliminator.cpp => ConditionalStateEliminator.cpp} (50%) rename src/solver/stateelimination/{ConditionalEliminator.h => ConditionalStateEliminator.h} (68%) create mode 100644 test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.cpp create mode 100644 test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.hpp create mode 100644 test/functional/solver/EliminationLinearEquationSolverTest.cpp diff --git a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseCore/SparseMatrix.h b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseCore/SparseMatrix.h index ba5e3a9b6..ddd0fc08c 100644 --- a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseCore/SparseMatrix.h +++ b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseCore/SparseMatrix.h @@ -384,7 +384,7 @@ class SparseMatrix eigen_assert( (m_outerIndex[outer+1]-m_outerIndex[outer]==0 || m_data.index(m_data.size()-1) Base; public: - SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) + SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1),m_detPermR(1) { initperfvalues(); } diff --git a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_column_bmod.h b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_column_bmod.h index cacc7e987..3731b299c 100644 --- a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_column_bmod.h +++ b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_column_bmod.h @@ -128,7 +128,7 @@ Index SparseLUImpl::column_bmod(const Index jcol, const Index nseg { irow = glu.lsub(isub); glu.lusup(nextlu) = dense(irow); - dense(irow) = Scalar(0.0); + dense(irow) = Scalar(0); ++nextlu; } diff --git a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h index 170610d9f..92920fd92 100644 --- a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h +++ b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h @@ -86,7 +86,7 @@ Index SparseLUImpl::copy_to_ucol(const Index jcol, const Index nse irow = glu.lsub(isub); glu.usub(nextu) = perm_r(irow); // Unlike the L part, the U part is stored in its final order glu.ucol(nextu) = dense(irow); - dense(irow) = Scalar(0.0); + dense(irow) = Scalar(0); nextu++; isub++; } diff --git a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_panel_bmod.h index 9d2ff2906..c79d9ff0c 100644 --- a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_panel_bmod.h +++ b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_panel_bmod.h @@ -122,7 +122,7 @@ void SparseLUImpl::panel_bmod(const Index m, const Index w, const Index isub = lptr + no_zeros; Index off = u_rows-segsize; - for (Index i = 0; i < off; i++) U(i,u_col) = 0; + for (Index i = 0; i < off; i++) U(i,u_col) = Scalar(0); for (Index i = 0; i < segsize; i++) { Index irow = glu.lsub(isub); @@ -172,7 +172,7 @@ void SparseLUImpl::panel_bmod(const Index m, const Index w, const { Index irow = glu.lsub(isub++); dense_col(irow) = U.coeff(i+off,u_col); - U.coeffRef(i+off,u_col) = 0; + U.coeffRef(i+off,u_col) = Scalar(0); } // Scatter l into SPA dense[] @@ -180,7 +180,7 @@ void SparseLUImpl::panel_bmod(const Index m, const Index w, const { Index irow = glu.lsub(isub++); dense_col(irow) -= L.coeff(i,u_col); - L.coeffRef(i,u_col) = 0; + L.coeffRef(i,u_col) = Scalar(0); } u_col++; } diff --git a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_pivotL.h b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_pivotL.h index 2e49ef667..4e4a2463e 100644 --- a/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_pivotL.h +++ b/resources/3rdparty/eigen-3.2.6/Eigen/src/SparseLU/SparseLU_pivotL.h @@ -71,42 +71,43 @@ Index SparseLUImpl::pivotL(const Index jcol, const RealScalar& dia // Determine the largest abs numerical value for partial pivoting Index diagind = iperm_c(jcol); // diagonal index - RealScalar pivmax(-1.0); - Index pivptr = nsupc; +// RealScalar pivmax(-1); + Index pivptr = nsupc; Index diag = emptyIdxLU; - RealScalar rtemp; - Index isub, icol, itemp, k; +// RealScalar rtemp; + Index isub, icol, itemp, k; for (isub = nsupc; isub < nsupr; ++isub) { - using std::abs; - rtemp = abs(lu_col_ptr[isub]); - if (rtemp > pivmax) { - pivmax = rtemp; +// using std::abs; +// rtemp = abs(lu_col_ptr[isub]); +// if (rtemp > pivmax) { +// pivmax = rtemp; pivptr = isub; - } +// } if (lsub_ptr[isub] == diagind) diag = isub; } // Test for singularity - if ( pivmax <= RealScalar(0.0) ) { - // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero - pivrow = pivmax < RealScalar(0.0) ? diagind : lsub_ptr[pivptr]; - perm_r(pivrow) = jcol; - return (jcol+1); - } +// if ( pivmax <= RealScalar(0) ) { +// // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero +// pivrow = pivmax < RealScalar(0) ? diagind : lsub_ptr[pivptr]; +// perm_r(pivrow) = jcol; +// return (jcol+1); +// } - RealScalar thresh = diagpivotthresh * pivmax; +// RealScalar thresh = diagpivotthresh * pivmax; // Choose appropriate pivotal element { // Test if the diagonal element can be used as a pivot (given the threshold value) - if (diag >= 0 ) - { +// if (diag >= 0 ) +// { // Diagonal element exists - using std::abs; - rtemp = abs(lu_col_ptr[diag]); - if (rtemp != 0.0 && rtemp >= thresh) pivptr = diag; - } +// using std::abs; +// rtemp = abs(lu_col_ptr[diag]); +// if (rtemp != 0.0 && rtemp >= thresh) + pivptr = diag; +// } pivrow = lsub_ptr[pivptr]; } @@ -125,7 +126,7 @@ Index SparseLUImpl::pivotL(const Index jcol, const RealScalar& dia } } // cdiv operations - Scalar temp = Scalar(1.0) / lu_col_ptr[nsupc]; + Scalar temp = Scalar(1) / lu_col_ptr[nsupc]; for (k = nsupc+1; k < nsupr; k++) lu_col_ptr[k] *= temp; return 0; diff --git a/src/adapters/EigenAdapter.cpp b/src/adapters/EigenAdapter.cpp index 45f5206bd..15be8c6ac 100644 --- a/src/adapters/EigenAdapter.cpp +++ b/src/adapters/EigenAdapter.cpp @@ -21,5 +21,7 @@ namespace storm { } template std::unique_ptr> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix const& matrix); + template std::unique_ptr> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix const& matrix); + template std::unique_ptr> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix const& matrix); } } \ No newline at end of file diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp index f6cb7e9c4..02cf81343 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp @@ -19,7 +19,7 @@ #include "src/logic/FragmentSpecification.h" #include "src/solver/stateelimination/LongRunAverageEliminator.h" -#include "src/solver/stateelimination/ConditionalEliminator.h" +#include "src/solver/stateelimination/ConditionalStateEliminator.h" #include "src/solver/stateelimination/PrioritizedStateEliminator.h" #include "src/solver/stateelimination/StaticStatePriorityQueue.h" #include "src/solver/stateelimination/DynamicStatePriorityQueue.h" @@ -685,7 +685,7 @@ namespace storm { STORM_LOG_INFO("Eliminating " << numberOfStatesToEliminate << " states using the state elimination technique." << std::endl); performPrioritizedStateElimination(statePriorities, flexibleMatrix, flexibleBackwardTransitions, oneStepProbabilities, this->getModel().getInitialStates(), true); - storm::solver::stateelimination::ConditionalEliminator stateEliminator = storm::solver::stateelimination::ConditionalEliminator(flexibleMatrix, flexibleBackwardTransitions, oneStepProbabilities, phiStates, psiStates); + storm::solver::stateelimination::ConditionalStateEliminator stateEliminator = storm::solver::stateelimination::ConditionalStateEliminator(flexibleMatrix, flexibleBackwardTransitions, oneStepProbabilities, phiStates, psiStates); // Eliminate the transitions going into the initial state (if there are any). if (!flexibleBackwardTransitions.getRow(*newInitialStates.begin()).empty()) { @@ -948,7 +948,7 @@ namespace storm { // Finally, eliminate the entry states (if we are required to do so). if (eliminateEntryStates) { STORM_LOG_TRACE("Finally, eliminating entry states."); - std::shared_ptr naivePriorities = createNaivePriorityQueue(entryStates); + std::shared_ptr naivePriorities = createStatePriorityQueue(entryStates); performPrioritizedStateElimination(naivePriorities, matrix, backwardTransitions, values, initialStates, computeResultsForInitialStatesOnly); STORM_LOG_TRACE("Eliminated/added entry states."); } else { diff --git a/src/settings/modules/MarkovChainSettings.cpp b/src/settings/modules/MarkovChainSettings.cpp index 589637103..2db5dbca7 100644 --- a/src/settings/modules/MarkovChainSettings.cpp +++ b/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 linearEquationSolver = {"gmm++", "native", "eigen"}; + std::vector linearEquationSolver = {"gmm++", "native", "eigen", "elimination"}; 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++, native and eigen.").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, eigen, elimination.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(linearEquationSolver)).setDefaultValueString("gmm++").build()).build()); std::vector ddLibraries = {"cudd", "sylvan"}; this->addOption(storm::settings::OptionBuilder(moduleName, ddLibraryOptionName, false, "Sets which library is preferred for decision-diagram operations.") @@ -87,6 +87,8 @@ namespace storm { return storm::solver::EquationSolverType::Native; } else if (equationSolverName == "eigen") { return storm::solver::EquationSolverType::Eigen; + } else if (equationSolverName == "elimination") { + return storm::solver::EquationSolverType::Elimination; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown equation solver '" << equationSolverName << "'."); } diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp index 7a6c83231..e364f1e19 100644 --- a/src/solver/EigenLinearEquationSolver.cpp +++ b/src/solver/EigenLinearEquationSolver.cpp @@ -9,7 +9,7 @@ namespace storm { namespace solver { template - EigenLinearEquationSolver::EigenLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix(A)), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner) { + EigenLinearEquationSolver::EigenLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix(A)), method(method), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner) { // Intentionally left empty. } @@ -44,31 +44,31 @@ namespace storm { template void EigenLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { // Translate the vectors x and b into Eigen's format. - Eigen::VectorXd eigenB(b.size()); + Eigen::Matrix eigenB(b.size()); for (uint64_t index = 0; index < b.size(); ++index) { eigenB[index] = b[index]; } - Eigen::VectorXd eigenX(x.size()); + Eigen::Matrix eigenX(x.size()); for (uint64_t index = 0; index < x.size(); ++index) { eigenX[index] = x[index]; } if (method == SolutionMethod::SparseLU) { - Eigen::SparseLU, Eigen::COLAMDOrdering> solver; + Eigen::SparseLU, Eigen::COLAMDOrdering> solver; solver.compute(*eigenA); solver._solve(eigenB, eigenX); } else { if (preconditioner == Preconditioner::Ilu) { - Eigen::BiCGSTAB, Eigen::IncompleteLUT> solver; + Eigen::BiCGSTAB, Eigen::IncompleteLUT> solver; solver.compute(*eigenA); solver._solve(eigenB, eigenX); } else if (preconditioner == Preconditioner::Diagonal) { - Eigen::BiCGSTAB, Eigen::DiagonalPreconditioner> solver; + Eigen::BiCGSTAB, Eigen::DiagonalPreconditioner> solver; solver.compute(*eigenA); solver._solve(eigenB, eigenX); } else { - Eigen::BiCGSTAB, Eigen::IdentityPreconditioner> solver; + Eigen::BiCGSTAB, Eigen::IdentityPreconditioner> solver; solver.compute(*eigenA); solver._solve(eigenB, eigenX); } @@ -82,16 +82,73 @@ namespace storm { template void EigenLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { + // Translate the vectors x and b into Eigen's format. + std::unique_ptr> eigenB; + if (b != nullptr) { + eigenB = std::make_unique>(b->size()); + for (uint64_t index = 0; index < b->size(); ++index) { + (*eigenB)[index] = (*b)[index]; + } + } + Eigen::Matrix 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]; + } + } + + // Specialization form storm::RationalNumber + + EigenLinearEquationSolver::EigenLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix(A)), method(method) { + // Intentionally left empty. + } + + void EigenLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + // Translate the vectors x and b into Eigen's format. + Eigen::Matrix eigenB(b.size()); + for (uint64_t index = 0; index < b.size(); ++index) { + eigenB[index] = b[index]; + } + + Eigen::Matrix eigenX(x.size()); + for (uint64_t index = 0; index < x.size(); ++index) { + eigenX[index] = x[index]; + } + if (method == SolutionMethod::SparseLU) { + Eigen::SparseLU, Eigen::COLAMDOrdering> 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]; + } + } + + void EigenLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { // Translate the vectors x and b into Eigen's format. - std::unique_ptr eigenB; + std::unique_ptr> eigenB; if (b != nullptr) { - eigenB = std::make_unique(b->size()); + eigenB = std::make_unique>(b->size()); for (uint64_t index = 0; index < b->size(); ++index) { (*eigenB)[index] = (*b)[index]; } } - Eigen::VectorXd eigenX(x.size()); + Eigen::Matrix eigenX(x.size()); for (uint64_t index = 0; index < x.size(); ++index) { eigenX[index] = x[index]; } @@ -110,8 +167,67 @@ namespace storm { } } - template class EigenLinearEquationSolver; + // Specialization form storm::RationalFunction + EigenLinearEquationSolver::EigenLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix(A)), method(method) { + // Intentionally left empty. + } + + void EigenLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + // Translate the vectors x and b into Eigen's format. + Eigen::Matrix eigenB(b.size()); + for (uint64_t index = 0; index < b.size(); ++index) { + eigenB[index] = b[index]; + } + + Eigen::Matrix eigenX(x.size()); + for (uint64_t index = 0; index < x.size(); ++index) { + eigenX[index] = x[index]; + } + + if (method == SolutionMethod::SparseLU) { + Eigen::SparseLU, Eigen::COLAMDOrdering> 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]; + } + } + + void EigenLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { + // Translate the vectors x and b into Eigen's format. + std::unique_ptr> eigenB; + if (b != nullptr) { + eigenB = std::make_unique>(b->size()); + for (uint64_t index = 0; index < b->size(); ++index) { + (*eigenB)[index] = (*b)[index]; + } + } + Eigen::Matrix 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; + template class EigenLinearEquationSolver; +// template class EigenLinearEquationSolver; } } \ No newline at end of file diff --git a/src/solver/EigenLinearEquationSolver.h b/src/solver/EigenLinearEquationSolver.h index 91e84f8f0..924a2af74 100644 --- a/src/solver/EigenLinearEquationSolver.h +++ b/src/solver/EigenLinearEquationSolver.h @@ -34,7 +34,7 @@ namespace storm { storm::storage::SparseMatrix const* originalA; // The (eigen) matrix associated with this equation solver. - std::unique_ptr> eigenA; + std::unique_ptr> eigenA; // The method to use for solving linear equation systems. SolutionMethod method; @@ -49,5 +49,52 @@ namespace storm { Preconditioner preconditioner; }; + template<> + class EigenLinearEquationSolver : public LinearEquationSolver { + public: + enum class SolutionMethod { + SparseLU + }; + + EigenLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method = SolutionMethod::SparseLU); + + virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; + + virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; + + private: + // A pointer to the original sparse matrix given to this solver. + storm::storage::SparseMatrix const* originalA; + + // The (eigen) matrix associated with this equation solver. + std::unique_ptr> eigenA; + + // The method to use for solving linear equation systems. + SolutionMethod method; + }; + + template<> + class EigenLinearEquationSolver : public LinearEquationSolver { + public: + enum class SolutionMethod { + SparseLU + }; + + EigenLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method = SolutionMethod::SparseLU); + + virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; + + virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; + + private: + // A pointer to the original sparse matrix given to this solver. + storm::storage::SparseMatrix const* originalA; + + // The (eigen) matrix associated with this equation solver. + std::unique_ptr> eigenA; + + // The method to use for solving linear equation systems. + SolutionMethod method; + }; } } \ No newline at end of file diff --git a/src/solver/EliminationLinearEquationSolver.cpp b/src/solver/EliminationLinearEquationSolver.cpp index 67a281b70..72e675569 100644 --- a/src/solver/EliminationLinearEquationSolver.cpp +++ b/src/solver/EliminationLinearEquationSolver.cpp @@ -1,9 +1,19 @@ #include "src/solver/EliminationLinearEquationSolver.h" +#include + +#include "src/solver/stateelimination/StatePriorityQueue.h" +#include "src/solver/stateelimination/PrioritizedStateEliminator.h" + +#include "src/utility/macros.h" #include "src/utility/vector.h" +#include "src/utility/stateelimination.h" namespace storm { namespace solver { + + using namespace stateelimination; + template EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix const& A) : A(A) { // Intentionally left empty. @@ -11,7 +21,54 @@ namespace storm { template void EliminationLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { - // TODO: implement state-elimination here. + STORM_LOG_WARN_COND(multiplyResult == nullptr, "Providing scratch memory will not improve the performance of this solver."); + +// std::cout << "input:" << std::endl; +// std::cout << "A:" << std::endl; +// std::cout << A << std::endl; +// std::cout << "b:" << std::endl; +// for (auto const& e : b) { +// std::cout << e << std::endl; +// } + + // Create a naive priority queue. + // TODO: improve. + std::vector allRows(x.size()); + std::iota(allRows.begin(), allRows.end(), 0); + std::shared_ptr priorityQueue = storm::utility::stateelimination::createStatePriorityQueue(allRows); + + // Initialize the solution to the right-hand side of the equation system. + x = b; + + // Translate the matrix and its transpose into the flexible format. + // We need to revert the matrix from the equation system format to the probability matrix format. That is, + // from (I-P), we construct P. Note that for the backwards transitions, this does not need to be done, as all + // entries are set to one anyway. + storm::storage::FlexibleSparseMatrix flexibleMatrix(A, false, true); + storm::storage::FlexibleSparseMatrix flexibleBackwardTransitions(A.transpose(), true, true); + +// std::cout << "intermediate:" << std::endl; +// std::cout << "flexibleMatrix:" << std::endl; +// std::cout << flexibleMatrix << std::endl; +// std::cout << "backward:" << std::endl; +// std::cout << flexibleBackwardTransitions << std::endl; + + // Create a state eliminator to perform the actual elimination. + PrioritizedStateEliminator eliminator(flexibleMatrix, flexibleBackwardTransitions, priorityQueue, x); + + std::cout << "eliminating" << std::endl; + while (priorityQueue->hasNext()) { + auto state = priorityQueue->pop(); + eliminator.eliminateState(state, false); + } + +// std::cout << "output:" << std::endl; +// std::cout << "x:" << std::endl; +// for (auto const& e : x) { +// std::cout << e << std::endl; +// } + + std::cout << "done" << std::endl; } template @@ -52,8 +109,6 @@ namespace storm { } template class EliminationLinearEquationSolver; - - // TODO: make this work with the proper implementation of solveEquationSystem. template class EliminationLinearEquationSolver; template class EliminationLinearEquationSolver; diff --git a/src/solver/SolverSelectionOptions.h b/src/solver/SolverSelectionOptions.h index 49f28d1de..0bc3cbf4a 100644 --- a/src/solver/SolverSelectionOptions.h +++ b/src/solver/SolverSelectionOptions.h @@ -9,7 +9,7 @@ namespace storm { ExtendEnumsWithSelectionField(MinMaxTechnique, PolicyIteration, ValueIteration) ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk) - ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Topological) + ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination, Topological) ExtendEnumsWithSelectionField(SmtSolverType, Z3, Mathsat) } } diff --git a/src/solver/stateelimination/ConditionalEliminator.cpp b/src/solver/stateelimination/ConditionalStateEliminator.cpp similarity index 50% rename from src/solver/stateelimination/ConditionalEliminator.cpp rename to src/solver/stateelimination/ConditionalStateEliminator.cpp index c45a31cb7..f15b3342b 100644 --- a/src/solver/stateelimination/ConditionalEliminator.cpp +++ b/src/solver/stateelimination/ConditionalStateEliminator.cpp @@ -1,4 +1,4 @@ -#include "src/solver/stateelimination/ConditionalEliminator.h" +#include "src/solver/stateelimination/ConditionalStateEliminator.h" #include "src/utility/macros.h" #include "src/utility/constants.h" @@ -8,26 +8,21 @@ namespace storm { namespace stateelimination { template - ConditionalEliminator::ConditionalEliminator(storm::storage::FlexibleSparseMatrix& transitionMatrix, storm::storage::FlexibleSparseMatrix& backwardTransitions, std::vector& oneStepProbabilities, storm::storage::BitVector& phiStates, storm::storage::BitVector& psiStates) : StateEliminator(transitionMatrix, backwardTransitions), oneStepProbabilities(oneStepProbabilities), phiStates(phiStates), psiStates(psiStates), filterLabel(StateLabel::NONE) { + ConditionalStateEliminator::ConditionalStateEliminator(storm::storage::FlexibleSparseMatrix& transitionMatrix, storm::storage::FlexibleSparseMatrix& backwardTransitions, std::vector& oneStepProbabilities, storm::storage::BitVector& phiStates, storm::storage::BitVector& psiStates) : StateEliminator(transitionMatrix, backwardTransitions), oneStepProbabilities(oneStepProbabilities), phiStates(phiStates), psiStates(psiStates), filterLabel(StateLabel::NONE) { } template - void ConditionalEliminator::updateValue(storm::storage::sparse::state_type const& state, ValueType const& loopProbability) { + void ConditionalStateEliminator::updateValue(storm::storage::sparse::state_type const& state, ValueType const& loopProbability) { oneStepProbabilities[state] = storm::utility::simplify(loopProbability * oneStepProbabilities[state]); } template - void ConditionalEliminator::updatePredecessor(storm::storage::sparse::state_type const& predecessor, ValueType const& probability, storm::storage::sparse::state_type const& state) { + void ConditionalStateEliminator::updatePredecessor(storm::storage::sparse::state_type const& predecessor, ValueType const& probability, storm::storage::sparse::state_type const& state) { oneStepProbabilities[predecessor] = storm::utility::simplify(oneStepProbabilities[predecessor] * storm::utility::simplify(probability * oneStepProbabilities[state])); } - - template - void ConditionalEliminator::updatePriority(storm::storage::sparse::state_type const& state) { - // Do nothing - } - + template - bool ConditionalEliminator::filterPredecessor(storm::storage::sparse::state_type const& state) { + bool ConditionalStateEliminator::filterPredecessor(storm::storage::sparse::state_type const& state) { // TODO find better solution than flag switch (filterLabel) { case StateLabel::PHI: @@ -41,34 +36,34 @@ namespace storm { } template - bool ConditionalEliminator::isFilterPredecessor() const { + bool ConditionalStateEliminator::isFilterPredecessor() const { return true; } template - void ConditionalEliminator::setFilterPhi() { + void ConditionalStateEliminator::setFilterPhi() { filterLabel = StateLabel::PHI; } template - void ConditionalEliminator::setFilterPsi() { + void ConditionalStateEliminator::setFilterPsi() { filterLabel = StateLabel::PSI; } template - void ConditionalEliminator::setFilter(StateLabel const& stateLabel) { + void ConditionalStateEliminator::setFilter(StateLabel const& stateLabel) { filterLabel = stateLabel; } template - void ConditionalEliminator::unsetFilter() { + void ConditionalStateEliminator::unsetFilter() { filterLabel = StateLabel::NONE; } - template class ConditionalEliminator; + template class ConditionalStateEliminator; #ifdef STORM_HAVE_CARL - template class ConditionalEliminator; + template class ConditionalStateEliminator; #endif } // namespace stateelimination diff --git a/src/solver/stateelimination/ConditionalEliminator.h b/src/solver/stateelimination/ConditionalStateEliminator.h similarity index 68% rename from src/solver/stateelimination/ConditionalEliminator.h rename to src/solver/stateelimination/ConditionalStateEliminator.h index ead6afeec..0d7cdcd5b 100644 --- a/src/solver/stateelimination/ConditionalEliminator.h +++ b/src/solver/stateelimination/ConditionalStateEliminator.h @@ -1,5 +1,5 @@ -#ifndef STORM_SOLVER_STATEELIMINATION_CONDITIONALELIMINATOR_H_ -#define STORM_SOLVER_STATEELIMINATION_CONDITIONALELIMINATOR_H_ +#ifndef STORM_SOLVER_STATEELIMINATION_CONDITIONALSTATEELIMINATOR_H_ +#define STORM_SOLVER_STATEELIMINATION_CONDITIONALSTATEELIMINATOR_H_ #include "src/solver/stateelimination/StateEliminator.h" @@ -10,16 +10,15 @@ namespace storm { namespace stateelimination { template - class ConditionalEliminator : public StateEliminator { + class ConditionalStateEliminator : public StateEliminator { public: enum class StateLabel { NONE, PHI, PSI }; - ConditionalEliminator(storm::storage::FlexibleSparseMatrix& transitionMatrix, storm::storage::FlexibleSparseMatrix& backwardTransitions, std::vector& oneStepProbabilities, storm::storage::BitVector& phiStates, storm::storage::BitVector& psiStates); + ConditionalStateEliminator(storm::storage::FlexibleSparseMatrix& transitionMatrix, storm::storage::FlexibleSparseMatrix& backwardTransitions, std::vector& oneStepProbabilities, storm::storage::BitVector& phiStates, storm::storage::BitVector& psiStates); // Instantiaton of Virtual methods void updateValue(storm::storage::sparse::state_type const& state, ValueType const& loopProbability) override; void updatePredecessor(storm::storage::sparse::state_type const& predecessor, ValueType const& probability, storm::storage::sparse::state_type const& state) override; - void updatePriority(storm::storage::sparse::state_type const& state) override; bool filterPredecessor(storm::storage::sparse::state_type const& state) override; bool isFilterPredecessor() const override; @@ -39,4 +38,4 @@ namespace storm { } // namespace storage } // namespace storm -#endif // STORM_SOLVER_STATEELIMINATION_CONDITIONALELIMINATOR_H_ +#endif // STORM_SOLVER_STATEELIMINATION_CONDITIONALSTATEELIMINATOR_H_ diff --git a/src/solver/stateelimination/PrioritizedStateEliminator.cpp b/src/solver/stateelimination/PrioritizedStateEliminator.cpp index b5915b659..5fa1cfc8a 100644 --- a/src/solver/stateelimination/PrioritizedStateEliminator.cpp +++ b/src/solver/stateelimination/PrioritizedStateEliminator.cpp @@ -29,7 +29,8 @@ namespace storm { } template class PrioritizedStateEliminator; - + template class PrioritizedStateEliminator; + #ifdef STORM_HAVE_CARL template class PrioritizedStateEliminator; #endif diff --git a/src/solver/stateelimination/StateEliminator.cpp b/src/solver/stateelimination/StateEliminator.cpp index eebe90bbc..ef4bb355e 100644 --- a/src/solver/stateelimination/StateEliminator.cpp +++ b/src/solver/stateelimination/StateEliminator.cpp @@ -264,6 +264,7 @@ namespace storm { } template class StateEliminator; + template class StateEliminator; #ifdef STORM_HAVE_CARL template class StateEliminator; diff --git a/src/storage/FlexibleSparseMatrix.cpp b/src/storage/FlexibleSparseMatrix.cpp index 2e3a4e876..563f86633 100644 --- a/src/storage/FlexibleSparseMatrix.cpp +++ b/src/storage/FlexibleSparseMatrix.cpp @@ -3,7 +3,10 @@ #include "src/storage/SparseMatrix.h" #include "src/storage/BitVector.h" #include "src/adapters/CarlAdapter.h" + +#include "src/utility/macros.h" #include "src/utility/constants.h" +#include "src/exceptions/InvalidArgumentException.h" namespace storm { namespace storage { @@ -13,7 +16,9 @@ namespace storm { } template - FlexibleSparseMatrix::FlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne) : data(matrix.getRowCount()), columnCount(matrix.getColumnCount()), nonzeroEntryCount(matrix.getNonzeroEntryCount()), trivialRowGrouping(matrix.hasTrivialRowGrouping()) { + FlexibleSparseMatrix::FlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne, bool revertEquationSystem) : data(matrix.getRowCount()), columnCount(matrix.getColumnCount()), nonzeroEntryCount(matrix.getNonzeroEntryCount()), trivialRowGrouping(matrix.hasTrivialRowGrouping()) { + STORM_LOG_THROW(!revertEquationSystem || trivialRowGrouping, storm::exceptions::InvalidArgumentException, "Illegal option for creating flexible matrix."); + if (!trivialRowGrouping) { rowGroupIndices = matrix.getRowGroupIndices(); } @@ -23,12 +28,31 @@ namespace storm { for (auto const& element : row) { // If the probability is zero, we skip this entry. if (storm::utility::isZero(element.getValue())) { - continue; + if (revertEquationSystem && rowIndex == element.getColumn()) { + getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one()); + } else { + continue; + } } if (setAllValuesToOne) { - getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one()); + if (revertEquationSystem && element.getColumn() == rowIndex && storm::utility::isOne(element.getValue())) { + continue; + } else { + getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one()); + } } else { - getRow(rowIndex).emplace_back(element); + if (revertEquationSystem) { + if (element.getColumn() == rowIndex) { + if (storm::utility::isOne(element.getValue())) { + continue; + } + getRow(rowIndex).emplace_back(element.getColumn(), storm::utility::one() - element.getValue()); + } else { + getRow(rowIndex).emplace_back(element.getColumn(), -element.getValue()); + } + } else { + getRow(rowIndex).emplace_back(element); + } } } } @@ -245,7 +269,9 @@ namespace storm { // Explicitly instantiate the matrix. template class FlexibleSparseMatrix; template std::ostream& operator<<(std::ostream& out, FlexibleSparseMatrix const& matrix); - + template class FlexibleSparseMatrix; + template std::ostream& operator<<(std::ostream& out, FlexibleSparseMatrix const& matrix); + #ifdef STORM_HAVE_CARL template class FlexibleSparseMatrix; template std::ostream& operator<<(std::ostream& out, FlexibleSparseMatrix const& matrix); diff --git a/src/storage/FlexibleSparseMatrix.h b/src/storage/FlexibleSparseMatrix.h index d5e6c2124..e2a1a5376 100644 --- a/src/storage/FlexibleSparseMatrix.h +++ b/src/storage/FlexibleSparseMatrix.h @@ -44,8 +44,10 @@ namespace storm { * Constructs a flexible sparse matrix from a sparse matrix. * @param matrix Sparse matrix to construct from. * @param setAllValuesToOne If true, all set entries are set to one. Default is false. + * @param revertEquationSystem If true, the matrix that will be created is the matrix (1-A), where A is the + * provided matrix. */ - FlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne = false); + FlexibleSparseMatrix(storm::storage::SparseMatrix const& matrix, bool setAllValuesToOne = false, bool revertEquationSystem = false); /*! * Reserves space for elements in row. diff --git a/src/utility/constants.cpp b/src/utility/constants.cpp index c3184c194..b2b182247 100644 --- a/src/utility/constants.cpp +++ b/src/utility/constants.cpp @@ -231,6 +231,24 @@ namespace storm { template storm::storage::MatrixEntry& simplify(storm::storage::MatrixEntry& matrixEntry); template storm::storage::MatrixEntry&& simplify(storm::storage::MatrixEntry&& matrixEntry); + // Instantiations for rational number. + template bool isOne(storm::RationalNumber const& value); + template bool isZero(storm::RationalNumber const& value); + template bool isConstant(storm::RationalNumber const& value); + + template storm::RationalNumber one(); + template storm::RationalNumber zero(); + + template double convertNumber(storm::RationalNumber const& number); + template storm::RationalNumber convertNumber(double const& number); + +// template storm::RationalNumber pow(storm::RationalNumber const& value, uint_fast64_t exponent); + + template storm::RationalNumber simplify(storm::RationalNumber value); + template storm::storage::MatrixEntry simplify(storm::storage::MatrixEntry matrixEntry); + template storm::storage::MatrixEntry& simplify(storm::storage::MatrixEntry& matrixEntry); + template storm::storage::MatrixEntry&& simplify(storm::storage::MatrixEntry&& matrixEntry); + #ifdef STORM_HAVE_CARL template bool isOne(RationalFunction const& value); template bool isZero(RationalFunction const& value); @@ -248,16 +266,6 @@ namespace storm { template Polynomial one(); template Polynomial zero(); - template bool isOne(RationalNumber const& value); - template bool isZero(RationalNumber const& value); - template bool isConstant(RationalNumber const& value); - - template RationalNumber one(); - template RationalNumber zero(); - - template double convertNumber(RationalNumber const& number); - template RationalNumber convertNumber(double const& number); - template bool isOne(Interval const& value); template bool isZero(Interval const& value); template bool isConstant(Interval const& value); diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp index 6c1d56077..58ef940f7 100644 --- a/src/utility/solver.cpp +++ b/src/utility/solver.cpp @@ -8,6 +8,7 @@ #include "src/solver/NativeLinearEquationSolver.h" #include "src/solver/GmmxxLinearEquationSolver.h" #include "src/solver/EigenLinearEquationSolver.h" +#include "src/solver/EliminationLinearEquationSolver.h" #include "src/solver/GameSolver.h" #include "src/solver/NativeMinMaxLinearEquationSolver.h" @@ -49,13 +50,18 @@ namespace storm { std::unique_ptr> LinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { storm::solver::EquationSolverType equationSolver = storm::settings::getModule().getEquationSolver(); switch (equationSolver) { - case storm::solver::EquationSolverType::Gmmxx: return std::unique_ptr>(new storm::solver::GmmxxLinearEquationSolver(matrix)); - case storm::solver::EquationSolverType::Native: return std::unique_ptr>(new storm::solver::NativeLinearEquationSolver(matrix)); - case storm::solver::EquationSolverType::Eigen: return std::unique_ptr>(new storm::solver::EigenLinearEquationSolver(matrix)); - default: return std::unique_ptr>(new storm::solver::GmmxxLinearEquationSolver(matrix)); + case storm::solver::EquationSolverType::Gmmxx: return std::make_unique>(matrix); + case storm::solver::EquationSolverType::Native: return std::make_unique>(matrix); + case storm::solver::EquationSolverType::Eigen: return std::make_unique>(matrix); + case storm::solver::EquationSolverType::Elimination: return std::make_unique>(matrix); + default: return std::make_unique>(matrix); } } + std::unique_ptr> LinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { + return std::make_unique>(matrix); + } + template std::unique_ptr> GmmxxLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { return std::unique_ptr>(new storm::solver::GmmxxLinearEquationSolver(matrix)); @@ -205,6 +211,7 @@ namespace storm { template class SymbolicGameSolverFactory; template class LinearEquationSolverFactory; template class GmmxxLinearEquationSolverFactory; + template class EigenLinearEquationSolverFactory; template class NativeLinearEquationSolverFactory; template class MinMaxLinearEquationSolverFactory; template class GameSolverFactory; diff --git a/src/utility/solver.h b/src/utility/solver.h index c0770605f..201cfccb1 100644 --- a/src/utility/solver.h +++ b/src/utility/solver.h @@ -4,8 +4,10 @@ #include #include #include -#include +#include "src/adapters/CarlAdapter.h" + +#include "src/storage/sparse/StateType.h" #include "src/storage/dd/DdType.h" #include "src/solver/SolverSelectionOptions.h" @@ -87,6 +89,18 @@ namespace storm { virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const; }; + template<> + class LinearEquationSolverFactory { + public: + /*! + * Creates a new linear equation solver instance with the given matrix. + * + * @param matrix The matrix that defines the equation system. + * @return A pointer to the newly created solver. + */ + virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const; + }; + template class NativeLinearEquationSolverFactory : public LinearEquationSolverFactory { public: diff --git a/src/utility/stateelimination.cpp b/src/utility/stateelimination.cpp index 22d830e53..bfb967da2 100644 --- a/src/utility/stateelimination.cpp +++ b/src/utility/stateelimination.cpp @@ -138,16 +138,21 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Illegal elimination order selected."); } - std::shared_ptr createNaivePriorityQueue(storm::storage::BitVector const& states) { + std::shared_ptr createStatePriorityQueue(storm::storage::BitVector const& states) { std::vector sortedStates(states.begin(), states.end()); return std::make_shared(sortedStates); } + std::shared_ptr createStatePriorityQueue(std::vector const& states) { + return std::make_shared(states); + } + template uint_fast64_t estimateComplexity(double const& value); template std::shared_ptr createStatePriorityQueue(boost::optional> const& distanceBasedStatePriorities, storm::storage::FlexibleSparseMatrix const& transitionMatrix, storm::storage::FlexibleSparseMatrix const& backwardTransitions, std::vector& oneStepProbabilities, storm::storage::BitVector const& states); template uint_fast64_t computeStatePenalty(storm::storage::sparse::state_type const& state, storm::storage::FlexibleSparseMatrix const& transitionMatrix, storm::storage::FlexibleSparseMatrix const& backwardTransitions, std::vector const& oneStepProbabilities); template uint_fast64_t computeStatePenaltyRegularExpression(storm::storage::sparse::state_type const& state, storm::storage::FlexibleSparseMatrix const& transitionMatrix, storm::storage::FlexibleSparseMatrix const& backwardTransitions, std::vector const& oneStepProbabilities); + template uint_fast64_t estimateComplexity(storm::RationalNumber const& value); #ifdef STORM_HAVE_CARL template uint_fast64_t estimateComplexity(storm::RationalFunction const& value); diff --git a/src/utility/stateelimination.h b/src/utility/stateelimination.h index 8703c4b4d..13c8904eb 100644 --- a/src/utility/stateelimination.h +++ b/src/utility/stateelimination.h @@ -53,7 +53,8 @@ namespace storm { template std::shared_ptr createStatePriorityQueue(boost::optional> const& stateDistances, storm::storage::FlexibleSparseMatrix const& transitionMatrix, storm::storage::FlexibleSparseMatrix const& backwardTransitions, std::vector& oneStepProbabilities, storm::storage::BitVector const& states); - std::shared_ptr createNaivePriorityQueue(storm::storage::BitVector const& states); + std::shared_ptr createStatePriorityQueue(storm::storage::BitVector const& states); + std::shared_ptr createStatePriorityQueue(std::vector const& states); } } diff --git a/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.cpp new file mode 100644 index 000000000..0a339713b --- /dev/null +++ b/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.cpp @@ -0,0 +1,321 @@ +#include "gtest/gtest.h" +#include "storm-config.h" + +#include "src/parser/FormulaParser.h" +#include "src/logic/Formulas.h" +#include "src/utility/solver.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" +#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "src/settings/SettingsManager.h" +#include "src/settings/modules/GeneralSettings.h" +#include "src/settings/modules/EigenEquationSolverSettings.h" +#include "src/settings/modules/NativeEquationSolverSettings.h" +#include "src/settings/SettingMemento.h" +#include "src/parser/AutoParser.h" +#include "src/parser/PrismParser.h" +#include "src/builder/ExplicitPrismModelBuilder.h" + +TEST(EigenDtmcPrctlModelCheckerTest, Die) { + std::shared_ptr> abstractModel = storm::parser::AutoParser<>::parseModel(STORM_CPP_BASE_PATH "/examples/dtmc/die/die.tra", STORM_CPP_BASE_PATH "/examples/dtmc/die/die.lab", "", STORM_CPP_BASE_PATH "/examples/dtmc/die/die.coin_flips.trans.rew"); + + // A parser that we use for conveniently constructing the formulas. + storm::parser::FormulaParser formulaParser; + + ASSERT_EQ(abstractModel->getType(), storm::models::ModelType::Dtmc); + + std::shared_ptr> dtmc = abstractModel->as>(); + + ASSERT_EQ(dtmc->getNumberOfStates(), 13ull); + ASSERT_EQ(dtmc->getNumberOfTransitions(), 20ull); + + storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::unique_ptr>(new storm::utility::solver::EigenLinearEquationSolverFactory())); + + std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]"); + + std::unique_ptr result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1.0 / 6.0, quantitativeResult1[0], storm::settings::getModule().getPrecision()); + + formula = formulaParser.parseSingleFormulaFromString("P=? [F \"two\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1.0 / 6.0, quantitativeResult2[0], storm::settings::getModule().getPrecision()); + + formula = formulaParser.parseSingleFormulaFromString("P=? [F \"three\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult3 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1.0 / 6.0, quantitativeResult3[0], storm::settings::getModule().getPrecision()); + + formula = formulaParser.parseSingleFormulaFromString("R=? [F \"done\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult4 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(11.0 / 3.0, quantitativeResult4[0], storm::settings::getModule().getPrecision()); +} + +TEST(EigenDtmcPrctlModelCheckerTest, Crowds) { + std::shared_ptr> abstractModel = storm::parser::AutoParser<>::parseModel(STORM_CPP_BASE_PATH "/examples/dtmc/crowds/crowds5_5.tra", STORM_CPP_BASE_PATH "/examples/dtmc/crowds/crowds5_5.lab", "", ""); + + ASSERT_EQ(abstractModel->getType(), storm::models::ModelType::Dtmc); + + // A parser that we use for conveniently constructing the formulas. + storm::parser::FormulaParser formulaParser; + + std::shared_ptr> dtmc = abstractModel->as>(); + + ASSERT_EQ(8607ull, dtmc->getNumberOfStates()); + ASSERT_EQ(15113ull, dtmc->getNumberOfTransitions()); + + storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::unique_ptr>(new storm::utility::solver::EigenLinearEquationSolverFactory())); + + std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observe0Greater1\"]"); + + std::unique_ptr result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.3328800375801578281, quantitativeResult1[0], storm::settings::getModule().getPrecision()); + + formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observeIGreater1\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.1522194965, quantitativeResult2[0], storm::settings::getModule().getPrecision()); + + formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observeOnlyTrueSender\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult3 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.32153724292835045, quantitativeResult3[0], storm::settings::getModule().getPrecision()); +} + +TEST(EigenDtmcPrctlModelCheckerTest, SynchronousLeader) { + std::shared_ptr> abstractModel = storm::parser::AutoParser<>::parseModel(STORM_CPP_BASE_PATH "/examples/dtmc/synchronous_leader/leader4_8.tra", STORM_CPP_BASE_PATH "/examples/dtmc/synchronous_leader/leader4_8.lab", "", STORM_CPP_BASE_PATH "/examples/dtmc/synchronous_leader/leader4_8.pick.trans.rew"); + + ASSERT_EQ(abstractModel->getType(), storm::models::ModelType::Dtmc); + + // A parser that we use for conveniently constructing the formulas. + storm::parser::FormulaParser formulaParser; + + std::shared_ptr> dtmc = abstractModel->as>(); + + ASSERT_EQ(12400ull, dtmc->getNumberOfStates()); + ASSERT_EQ(16495ull, dtmc->getNumberOfTransitions()); + + storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::unique_ptr>(new storm::utility::solver::EigenLinearEquationSolverFactory())); + + std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("P=? [F \"elected\"]"); + + std::unique_ptr result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1.0, quantitativeResult1[0], storm::settings::getModule().getPrecision()); + + formula = formulaParser.parseSingleFormulaFromString("P=? [F<=20 \"elected\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.9999965911265462636, quantitativeResult2[0], storm::settings::getModule().getPrecision()); + + formula = formulaParser.parseSingleFormulaFromString("R=? [F \"elected\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult3 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1.0448979591836789, quantitativeResult3[0], storm::settings::getModule().getPrecision()); +} + +TEST(EigenDtmcPrctlModelCheckerTest, LRASingleBscc) { + storm::storage::SparseMatrixBuilder matrixBuilder; + std::shared_ptr> dtmc; + + // A parser that we use for conveniently constructing the formulas. + storm::parser::FormulaParser formulaParser; + + { + matrixBuilder = storm::storage::SparseMatrixBuilder(2, 2, 2); + matrixBuilder.addNextValue(0, 1, 1.); + matrixBuilder.addNextValue(1, 0, 1.); + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(2); + ap.addLabel("a"); + ap.addLabelToState("a", 1); + + dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap)); + + storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + + std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]"); + + std::unique_ptr result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::getModule().getPrecision()); + } + { + matrixBuilder = storm::storage::SparseMatrixBuilder(2, 2, 4); + matrixBuilder.addNextValue(0, 0, .5); + matrixBuilder.addNextValue(0, 1, .5); + matrixBuilder.addNextValue(1, 0, .5); + matrixBuilder.addNextValue(1, 1, .5); + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(2); + ap.addLabel("a"); + ap.addLabelToState("a", 1); + + dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap)); + + storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + + std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]"); + + std::unique_ptr result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::getModule().getPrecision()); + } + + { + matrixBuilder = storm::storage::SparseMatrixBuilder(3, 3, 3); + matrixBuilder.addNextValue(0, 1, 1); + matrixBuilder.addNextValue(1, 2, 1); + matrixBuilder.addNextValue(2, 0, 1); + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(3); + ap.addLabel("a"); + ap.addLabelToState("a", 2); + + dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap)); + + storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::unique_ptr>(new storm::utility::solver::EigenLinearEquationSolverFactory())); + + std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]"); + + std::unique_ptr result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::getModule().getPrecision()); + } +} + +TEST(EigenDtmcPrctlModelCheckerTest, LRA) { + storm::storage::SparseMatrixBuilder matrixBuilder; + std::shared_ptr> dtmc; + + // A parser that we use for conveniently constructing the formulas. + storm::parser::FormulaParser formulaParser; + + { + matrixBuilder = storm::storage::SparseMatrixBuilder(15, 15, 20, true); + matrixBuilder.addNextValue(0, 1, 1); + matrixBuilder.addNextValue(1, 4, 0.7); + matrixBuilder.addNextValue(1, 6, 0.3); + matrixBuilder.addNextValue(2, 0, 1); + + matrixBuilder.addNextValue(3, 5, 0.8); + matrixBuilder.addNextValue(3, 9, 0.2); + matrixBuilder.addNextValue(4, 3, 1); + matrixBuilder.addNextValue(5, 3, 1); + + matrixBuilder.addNextValue(6, 7, 1); + matrixBuilder.addNextValue(7, 8, 1); + matrixBuilder.addNextValue(8, 6, 1); + + matrixBuilder.addNextValue(9, 10, 1); + matrixBuilder.addNextValue(10, 9, 1); + matrixBuilder.addNextValue(11, 9, 1); + + matrixBuilder.addNextValue(12, 5, 0.4); + matrixBuilder.addNextValue(12, 8, 0.3); + matrixBuilder.addNextValue(12, 11, 0.3); + + matrixBuilder.addNextValue(13, 7, 0.7); + matrixBuilder.addNextValue(13, 12, 0.3); + + matrixBuilder.addNextValue(14, 12, 1); + + storm::storage::SparseMatrix transitionMatrix = matrixBuilder.build(); + + storm::models::sparse::StateLabeling ap(15); + ap.addLabel("a"); + ap.addLabelToState("a", 1); + ap.addLabelToState("a", 4); + ap.addLabelToState("a", 5); + ap.addLabelToState("a", 7); + ap.addLabelToState("a", 11); + ap.addLabelToState("a", 13); + ap.addLabelToState("a", 14); + + dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap)); + + storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::unique_ptr>(new storm::utility::solver::EigenLinearEquationSolverFactory())); + + std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]"); + + std::unique_ptr result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + + EXPECT_NEAR(0.3 / 3., quantitativeResult1[0], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult1[3], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(1. / 3., quantitativeResult1[6], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(0.0, quantitativeResult1[9], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(0.3 / 3., quantitativeResult1[12], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(.79 / 3., quantitativeResult1[13], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(0.3 / 3., quantitativeResult1[14], storm::settings::getModule().getPrecision()); + } +} + +TEST(EigenDtmcPrctlModelCheckerTest, Conditional) { + storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/modelchecker/test_conditional.pm"); + + std::shared_ptr> model = storm::builder::ExplicitPrismModelBuilder(program).translate(); + ASSERT_TRUE(model->getType() == storm::models::ModelType::Dtmc); + ASSERT_EQ(4ul, model->getNumberOfStates()); + ASSERT_EQ(5ul, model->getNumberOfTransitions()); + + std::shared_ptr> dtmc = model->as>(); + + storm::modelchecker::SparseDtmcPrctlModelChecker> checker(*dtmc, std::unique_ptr>(new storm::utility::solver::EigenLinearEquationSolverFactory())); + + // A parser that we use for conveniently constructing the formulas. + storm::parser::FormulaParser formulaParser; + std::shared_ptr formula = formulaParser.parseSingleFormulaFromString("P=? [F \"target\"]"); + + std::unique_ptr result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult1 = result->asExplicitQuantitativeCheckResult(); + EXPECT_NEAR(0.5, quantitativeResult1[0], storm::settings::getModule().getPrecision()); + + formula = formulaParser.parseSingleFormulaFromString("P=? [F \"target\" || F \"condition\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult2 = result->asExplicitQuantitativeCheckResult(); + EXPECT_NEAR(storm::utility::one(), quantitativeResult2[0], storm::settings::getModule().getPrecision()); + + formula = formulaParser.parseSingleFormulaFromString("R=? [F \"target\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult3 = result->asExplicitQuantitativeCheckResult(); + EXPECT_EQ(storm::utility::infinity(), quantitativeResult3[0]); + + formula = formulaParser.parseSingleFormulaFromString("R=? [F \"target\" || F \"condition\"]"); + + result = checker.check(*formula); + storm::modelchecker::ExplicitQuantitativeCheckResult& quantitativeResult4 = result->asExplicitQuantitativeCheckResult(); + EXPECT_NEAR(storm::utility::one(), quantitativeResult4[0], storm::settings::getModule().getPrecision()); +} diff --git a/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.hpp b/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.hpp new file mode 100644 index 000000000..ba9cac61b --- /dev/null +++ b/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.hpp @@ -0,0 +1,14 @@ +// +// EigenDtmcPrctlModelCheckerTest.hpp +// storm +// +// Created by Christian Dehnert on 19.06.16. +// +// + +#ifndef EigenDtmcPrctlModelCheckerTest_hpp +#define EigenDtmcPrctlModelCheckerTest_hpp + +#include + +#endif /* EigenDtmcPrctlModelCheckerTest_hpp */ diff --git a/test/functional/solver/EigenLinearEquationSolverTest.cpp b/test/functional/solver/EigenLinearEquationSolverTest.cpp index 87703de9a..d83032301 100644 --- a/test/functional/solver/EigenLinearEquationSolverTest.cpp +++ b/test/functional/solver/EigenLinearEquationSolverTest.cpp @@ -4,7 +4,8 @@ #include "src/solver/EigenLinearEquationSolver.h" #include "src/settings/SettingsManager.h" -#include "src/settings/modules/GmmxxEquationSolverSettings.h" +#include "src/utility/constants.h" +#include "src/settings/modules/EigenEquationSolverSettings.h" TEST(EigenLinearEquationSolver, SolveWithStandardOptions) { ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); @@ -29,9 +30,9 @@ TEST(EigenLinearEquationSolver, SolveWithStandardOptions) { storm::solver::EigenLinearEquationSolver solver(A); ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); - ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); } TEST(EigenLinearEquationSolver, SparseLU) { @@ -57,9 +58,65 @@ TEST(EigenLinearEquationSolver, SparseLU) { storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::SparseLU, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::None); ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); - ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); +} + +TEST(EigenLinearEquationSolver, SparseLU_RationalNumber) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder 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 A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(3); + std::vector b = {16, -4, -7}; + + ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::SparseLU)); + + storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::SparseLU); + ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); + ASSERT_TRUE(storm::utility::isOne(x[0])); + ASSERT_TRUE(x[1] == 3); + ASSERT_TRUE(x[2] == -1); +} + +TEST(EigenLinearEquationSolver, SparseLU_RationalFunction) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder builder; + ASSERT_NO_THROW(builder.addNextValue(0, 0, storm::RationalFunction(2))); + ASSERT_NO_THROW(builder.addNextValue(0, 1, storm::RationalFunction(4))); + ASSERT_NO_THROW(builder.addNextValue(0, 2, storm::RationalFunction(-2))); + ASSERT_NO_THROW(builder.addNextValue(1, 0, storm::RationalFunction(4))); + ASSERT_NO_THROW(builder.addNextValue(1, 1, storm::RationalFunction(-1))); + ASSERT_NO_THROW(builder.addNextValue(1, 2, storm::RationalFunction(5))); + ASSERT_NO_THROW(builder.addNextValue(2, 0, storm::RationalFunction(-1))); + ASSERT_NO_THROW(builder.addNextValue(2, 1, storm::RationalFunction(-1))); + ASSERT_NO_THROW(builder.addNextValue(2, 2, storm::RationalFunction(3))); + + storm::storage::SparseMatrix A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(3); + std::vector b = {storm::RationalFunction(16), storm::RationalFunction(-4), storm::RationalFunction(-7)}; + + ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::SparseLU)); + + storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::SparseLU); + ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); + ASSERT_TRUE(storm::utility::isOne(x[0])); + ASSERT_TRUE(x[1] == storm::RationalFunction(3)); + ASSERT_TRUE(x[2] == storm::RationalFunction(-1)); } TEST(EigenLinearEquationSolver, BiCGSTAB) { @@ -85,9 +142,9 @@ TEST(EigenLinearEquationSolver, BiCGSTAB) { storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::None); ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); - ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); } TEST(EigenLinearEquationSolver, BiCGSTAB_Ilu) { @@ -113,9 +170,9 @@ TEST(EigenLinearEquationSolver, BiCGSTAB_Ilu) { storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::Ilu); ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); - ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); } TEST(EigenLinearEquationSolver, BiCGSTAB_Diagonal) { @@ -141,9 +198,9 @@ TEST(EigenLinearEquationSolver, BiCGSTAB_Diagonal) { storm::solver::EigenLinearEquationSolver solver(A, storm::solver::EigenLinearEquationSolver::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver::Preconditioner::Diagonal); ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); - ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); - ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); } TEST(EigenLinearEquationSolver, MatrixVectorMultplication) { @@ -166,5 +223,5 @@ TEST(EigenLinearEquationSolver, MatrixVectorMultplication) { storm::solver::EigenLinearEquationSolver solver(A); ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(x, nullptr, 4)); - ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); } \ No newline at end of file diff --git a/test/functional/solver/EliminationLinearEquationSolverTest.cpp b/test/functional/solver/EliminationLinearEquationSolverTest.cpp new file mode 100644 index 000000000..39007a5d7 --- /dev/null +++ b/test/functional/solver/EliminationLinearEquationSolverTest.cpp @@ -0,0 +1,58 @@ +#include "gtest/gtest.h" +#include "storm-config.h" + +#include "src/solver/EliminationLinearEquationSolver.h" +#include "src/settings/SettingsManager.h" + +#include "src/settings/modules/GmmxxEquationSolverSettings.h" + +TEST(EliminationLinearEquationSolver, Solve) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder 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 A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(3); + std::vector b = {16, -4, -7}; + + ASSERT_NO_THROW(storm::solver::EliminationLinearEquationSolver solver(A)); + + storm::solver::EliminationLinearEquationSolver solver(A); + ASSERT_NO_THROW(solver.solveEquationSystem(x, b)); + ASSERT_LT(std::abs(x[0] - 1), 1e-15); + ASSERT_LT(std::abs(x[1] - 3), 1e-15); + ASSERT_LT(std::abs(x[2] - (-1)), 1e-15); +} + +TEST(EliminationLinearEquationSolver, MatrixVectorMultplication) { + ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); + storm::storage::SparseMatrixBuilder 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 A; + ASSERT_NO_THROW(A = builder.build()); + + std::vector x(5); + x[4] = 1; + + storm::solver::EliminationLinearEquationSolver solver(A); + ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(x, nullptr, 4)); + ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule().getPrecision()); +} \ No newline at end of file