Browse Source

added 'convergence' (rather success) checks for Eigen solver

Former-commit-id: 25a6fb1d77
tempestpy_adaptions
dehnert 9 years ago
parent
commit
b67e3d6e7b
  1. 18
      src/solver/EigenLinearEquationSolver.cpp
  2. 2
      src/solver/EigenLinearEquationSolver.h

18
src/solver/EigenLinearEquationSolver.cpp

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

2
src/solver/EigenLinearEquationSolver.h

@ -66,7 +66,7 @@ namespace storm {
virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override; virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override; virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
virtual bool solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override; virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
EigenLinearEquationSolverSettings<ValueType>& getSettings(); EigenLinearEquationSolverSettings<ValueType>& getSettings();

Loading…
Cancel
Save