From: Guido Kanschat Date: Thu, 29 Mar 2001 23:35:15 +0000 (+0000) Subject: solvers throw exception on failure X-Git-Tag: v8.0.0~19438 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e4197f96bff7b1ac69d6ab8770294a58f6d37656;p=dealii.git solvers throw exception on failure git-svn-id: https://svn.dealii.org/trunk@4342 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 41224c6513..917b44fa26 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -238,6 +238,14 @@ documentation, etc.

lac

    +
  1. Changed: Solver functions solve + are void now. If the solver has not converged within the + maximum number of iterations, it throws an exception of type + ExcNoConvergence. +
    + (GK 2001/03/29) +

    +
  2. New: The functions FullMatrix::mmult and * l2-norm. */ template - typename Solver::ReturnState + void solve (double &value, const MATRIX &A, VECTOR &x); @@ -218,7 +218,7 @@ class EigenInverse : private Solver * l2-norm. */ template - typename Solver::ReturnState + void solve (double &value, const MATRIX &A, VECTOR &x); @@ -301,7 +301,7 @@ EigenPower::~EigenPower () template template -typename Solver::ReturnState +void EigenPower::solve (double &value, const MATRIX &A, VECTOR &x) @@ -362,10 +362,9 @@ EigenPower::solve (double &value, deallog.pop(); // Output - if (conv == SolverControl::failure) - return exceeded; - else - return success; + AssertThrow(control().last_check() == SolverControl::success, + typename Solver::ExcNoConvergence(control().last_step(), + control().last_value())); } //----------------------------------------------------------------------// @@ -386,7 +385,7 @@ EigenInverse::~EigenInverse () template template -typename Solver::ReturnState +void EigenInverse::solve (double &value, const MATRIX &A, VECTOR &x) @@ -472,10 +471,9 @@ EigenInverse::solve (double &value, deallog.pop(); // Output - if (conv == SolverControl::failure) - return exceeded; - else - return success; + AssertThrow(control().last_check() == SolverControl::success, + typename Solver::ExcNoConvergence(control().last_step(), + control().last_value())); } #endif diff --git a/deal.II/lac/include/lac/solver.h b/deal.II/lac/include/lac/solver.h index f76a65bfff..daa52fa542 100644 --- a/deal.II/lac/include/lac/solver.h +++ b/deal.II/lac/include/lac/solver.h @@ -128,21 +128,12 @@ class SolverControl; * changes and it is still possible in a simple way to give these additional data to * the @p{SolverSelector} object for each solver which it may use. * - * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997, 1998, 1999 + * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997-2001 */ template > class Solver { public: - /** - * Declare possible return values of a - * solver object. - */ - enum ReturnState { - success=0, exceeded, breakdown - }; - - /** * Constructor. Assign a control * object which stores the required @@ -162,7 +153,11 @@ class Solver * Access to control object. */ SolverControl& control() const; - + + DeclException2(ExcNoConvergence, int, double, + << "Iteration did not converge after " << arg1 + << " steps. Final residual " << arg2); + protected: /** diff --git a/deal.II/lac/include/lac/solver_bicgstab.h b/deal.II/lac/include/lac/solver_bicgstab.h index 63a282e657..9e6676879c 100644 --- a/deal.II/lac/include/lac/solver_bicgstab.h +++ b/deal.II/lac/include/lac/solver_bicgstab.h @@ -97,7 +97,7 @@ class SolverBicgstab : public Subscriptor, private Solver * Solve primal problem only. */ template - typename Solver::ReturnState + void solve (const MATRIX &A, VECTOR &x, const VECTOR &b, @@ -212,7 +212,7 @@ class SolverBicgstab : public Subscriptor, private Solver * The iteration loop itself. */ template - typename Solver::ReturnState + bool iterate(const MATRIX& A, const PRECONDITIONER& precondition); }; @@ -280,7 +280,7 @@ SolverBicgstab::print_vectors(const unsigned int, template template -typename Solver::ReturnState +bool SolverBicgstab::iterate(const MATRIX& A, const PRECONDITIONER& precondition) { @@ -312,7 +312,7 @@ SolverBicgstab::iterate(const MATRIX& A, //TODO:[GK] Find better breakdown criterion if (fabs(alpha) > 1.e10) - return typename Solver::ReturnState(breakdown); + return true; s.equ(1., r, -alpha, v); precondition.vmult(z,s); @@ -322,18 +322,18 @@ SolverBicgstab::iterate(const MATRIX& A, Vx->add(alpha, y, omega, z); r.equ(1., s, -omega, t); - state = control().check(++step, criterion(A, *Vx, *Vb)); + res = criterion(A, *Vx, *Vb); + state = control().check(++step, res); print_vectors(step, *Vx, r, y); } while (state == SolverControl::iterate); - if (state == SolverControl::success) return success; - return exceeded; + return false; } template template -typename Solver::ReturnState +void SolverBicgstab::solve(const MATRIX &A, VECTOR &x, const VECTOR &b, @@ -354,7 +354,7 @@ SolverBicgstab::solve(const MATRIX &A, step = 0; - typename Solver::ReturnState state = breakdown; + bool state; do { @@ -363,7 +363,7 @@ SolverBicgstab::solve(const MATRIX &A, if (start(A) == SolverControl::success) break; state = iterate(A, precondition); } - while (state == breakdown); + while (state); memory.free(Vr); memory.free(Vrbar); @@ -375,7 +375,10 @@ SolverBicgstab::solve(const MATRIX &A, memory.free(Vv); deallog.pop(); - return state; + + AssertThrow(control().last_check() == SolverControl::success, + typename Solver::ExcNoConvergence(control().last_step(), + control().last_value())); } diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index 79c4369aca..bf81acefc4 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -87,7 +88,7 @@ class SolverCG : public Subscriptor, private Solver * Solver method. */ template - typename Solver::ReturnState + void solve (const MATRIX &A, VECTOR &x, const VECTOR &b, @@ -187,7 +188,7 @@ SolverCG::print_vectors(const unsigned int, template template -typename Solver::ReturnState +void SolverCG::solve (const MATRIX &A, VECTOR &x, const VECTOR &b, @@ -231,7 +232,7 @@ SolverCG::solve (const MATRIX &A, memory.free(Vz); memory.free(VAp); deallog.pop(); - return success; + return; }; g.scale(-1.); @@ -281,10 +282,9 @@ SolverCG::solve (const MATRIX &A, // Output deallog.pop(); - if (conv == SolverControl::failure) - return exceeded; - else - return success; + AssertThrow(control().last_check() == SolverControl::success, + typename Solver::ExcNoConvergence(control().last_step(), + control().last_value())); }; diff --git a/deal.II/lac/include/lac/solver_control.h b/deal.II/lac/include/lac/solver_control.h index 648d8fef70..839e11e2ff 100644 --- a/deal.II/lac/include/lac/solver_control.h +++ b/deal.II/lac/include/lac/solver_control.h @@ -153,6 +153,11 @@ class SolverControl : public Subscriptor virtual State check (const unsigned int step, const double check_value); + /** + * Return the result of the last check operation. + */ + State last_check() const; + /** * Return the convergence value of last * iteration step for which @p{check} was @@ -236,6 +241,11 @@ class SolverControl : public Subscriptor */ double tol; + /** + * Result of last check operation. + */ + State lcheck; + /** * Last value of the convergence criterion. */ diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index 44764f23ad..cef56c0066 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -110,10 +110,10 @@ class SolverGMRES : public Subscriptor, private Solver * Solver method. */ template - typename Solver::ReturnState solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER &precondition); + void solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER& precondition); DeclException1 (ExcTooFewTmpVectors, int, @@ -190,11 +190,11 @@ SolverGMRES::givens_rotation (Vector &h, template template -typename Solver::ReturnState -SolverGMRES::solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER &precondition) +void +SolverGMRES::solve (const MATRIX& A, + VECTOR & x, + const VECTOR& b, + const PRECONDITIONER& precondition) { // this code was written a very // long time ago by people not @@ -397,10 +397,10 @@ unsigned int dim = 0; memory.free (tmp_vectors[tmp]); deallog.pop(); - if (iteration_state) - return success; - else - return exceeded; + + AssertThrow(control().last_check() == SolverControl::success, + typename Solver::ExcNoConvergence(control().last_step(), + control().last_value())); }; diff --git a/deal.II/lac/include/lac/solver_minres.h b/deal.II/lac/include/lac/solver_minres.h index 163121e4b3..fa7f86742c 100644 --- a/deal.II/lac/include/lac/solver_minres.h +++ b/deal.II/lac/include/lac/solver_minres.h @@ -77,7 +77,7 @@ class SolverMinRes : public Subscriptor, private Solver * Solver method. */ template - typename Solver::ReturnState + void solve (const MATRIX &A, VECTOR &x, const VECTOR &b, @@ -351,10 +351,9 @@ SolverMinRes::solve (const MATRIX &A, // Output deallog.pop (); - if (conv == SolverControl::failure) - return exceeded; - - return success; + AssertThrow(control().last_check() == SolverControl::success, + typename Solver::ExcNoConvergence(control().last_step(), + control().last_value())); }; diff --git a/deal.II/lac/include/lac/solver_qmrs.h b/deal.II/lac/include/lac/solver_qmrs.h index fa6f5d1c6e..4aa2aedf9b 100644 --- a/deal.II/lac/include/lac/solver_qmrs.h +++ b/deal.II/lac/include/lac/solver_qmrs.h @@ -109,11 +109,11 @@ class SolverQMRS : public Subscriptor, private Solver * Solver method. */ template - typename Solver::ReturnState + void solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition); + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER& precondition); /** * Interface for derived class. @@ -175,7 +175,7 @@ class SolverQMRS : public Subscriptor, private Solver * The iteration loop itself. */ template - typename Solver::ReturnState + bool iterate(const MATRIX& A, const PRECONDITIONER& precondition); /** * The current iteration step. @@ -217,7 +217,7 @@ SolverQMRS::print_vectors(const unsigned int, template template -typename Solver::ReturnState +void SolverQMRS::solve (const MATRIX &A, VECTOR &x, const VECTOR &b, @@ -244,7 +244,7 @@ SolverQMRS::solve (const MATRIX &A, step = 0; - typename Solver::ReturnState state = breakdown; + bool state; do { @@ -252,7 +252,7 @@ SolverQMRS::solve (const MATRIX &A, deallog << "Restart step " << step << std::endl; state = iterate(A, precondition); } - while (state == breakdown); + while (state); // Deallocate Memory @@ -265,15 +265,17 @@ SolverQMRS::solve (const MATRIX &A, // Output deallog.pop(); - - return state; + + AssertThrow(control().last_check() == SolverControl::success, + typename Solver::ExcNoConvergence(control().last_step(), + control().last_value())); }; template template -typename Solver::ReturnState +bool SolverQMRS::iterate(const MATRIX& A, const PRECONDITIONER& precondition) { @@ -309,7 +311,7 @@ SolverQMRS::iterate(const MATRIX& A, res = v.l2_norm(); if (control().check(step, res) == SolverControl::success) - return success; + return false; p = v; @@ -331,7 +333,7 @@ while (state == SolverControl::iterate) //TODO:[GK] Find a really good breakdown criterion. The absolute one detects breakdown instead of convergence if (fabs(sigma/rho) < additional_data.breakdown) - return breakdown; + return true; // Step 3 alpha = rho/sigma; //deallog << "alpha:" << alpha << std::endl; @@ -361,13 +363,12 @@ while (state == SolverControl::iterate) else res = sqrt((it+1)*tau); state = control().check(step,res); - if (state == SolverControl::success) - return success; - else if (state == SolverControl::failure) - return exceeded; + if ((state == SolverControl::success) + || (state == SolverControl::failure)) + return false; // Step 6 if (fabs(rho) < additional_data.breakdown) - return breakdown; + return true; // Step 7 rho_old = rho; precondition.vmult(q,v); @@ -377,7 +378,7 @@ while (state == SolverControl::iterate) p.sadd(beta,v); precondition.vmult(q,p); } - return exceeded; + return false; } diff --git a/deal.II/lac/include/lac/solver_richardson.h b/deal.II/lac/include/lac/solver_richardson.h index 847f2b4097..103550ec44 100644 --- a/deal.II/lac/include/lac/solver_richardson.h +++ b/deal.II/lac/include/lac/solver_richardson.h @@ -79,19 +79,19 @@ class SolverRichardson : public Subscriptor, private Solver * Solve $Ax=b$ for $x$. */ template - typename Solver::ReturnState solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition); + void solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER& precondition); /** * Solve $A^Tx=b$ for $x$. */ template - typename Solver::ReturnState Tsolve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition); + void Tsolve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER& precondition); /** * Set the damping-coefficient. @@ -168,7 +168,7 @@ SolverRichardson::~SolverRichardson() template template -typename Solver::ReturnState +void SolverRichardson::solve (const MATRIX &A, VECTOR &x, const VECTOR &b, @@ -206,16 +206,15 @@ SolverRichardson::solve (const MATRIX &A, deallog.pop(); // Output - if (conv == SolverControl::failure) - return exceeded; - else - return success; + AssertThrow(control().last_check() == SolverControl::success, + typename Solver::ExcNoConvergence(control().last_step(), + control().last_value())); } template template -typename Solver::ReturnState +void SolverRichardson::Tsolve (const MATRIX &A, VECTOR &x, const VECTOR &b, @@ -253,10 +252,9 @@ SolverRichardson::Tsolve (const MATRIX &A, deallog.pop(); // Output - if (conv == SolverControl::failure) - return exceeded; - else - return success; + AssertThrow(control().last_check() == SolverControl::success, + typename Solver::ExcNoConvergence(control().last_step(), + control().last_value())); } diff --git a/deal.II/lac/include/lac/solver_selector.h b/deal.II/lac/include/lac/solver_selector.h index be805096f7..268c441e94 100644 --- a/deal.II/lac/include/lac/solver_selector.h +++ b/deal.II/lac/include/lac/solver_selector.h @@ -111,10 +111,10 @@ class SolverSelector : public Subscriptor * */ template - typename Solver::ReturnState solve(const Matrix &A, - Vector &x, - const Vector &b, - const Preconditioner &precond) const; + void solve(const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner &precond) const; /** * Set the additional data. For more @@ -220,7 +220,7 @@ SolverSelector::~SolverSelector() template template -typename Solver::ReturnState +void SolverSelector::solve(const Matrix &A, Vector &x, const Vector &b, @@ -230,30 +230,28 @@ SolverSelector::solve(const Matrix &A, { SolverRichardson solver(*control,*vector_memory, richardson_data); - return solver.solve(A,x,b,precond); + solver.solve(A,x,b,precond); } else if (solver_name=="cg") { SolverCG solver(*control,*vector_memory, cg_data); - return solver.solve(A,x,b,precond); + solver.solve(A,x,b,precond); } else if (solver_name=="bicgstab") { SolverBicgstab solver(*control,*vector_memory, bicgstab_data); - return solver.solve(A,x,b,precond); + solver.solve(A,x,b,precond); } else if (solver_name=="gmres") { SolverGMRES solver(*control,*vector_memory, gmres_data); - return solver.solve(A,x,b,precond); + solver.solve(A,x,b,precond); } else Assert(false,ExcSolverDoesNotExist(solver_name)); - - return Solver::breakdown; }; diff --git a/deal.II/lac/source/solver_control.cc b/deal.II/lac/source/solver_control.cc index e9c2e9c6c4..27f671efb7 100644 --- a/deal.II/lac/source/solver_control.cc +++ b/deal.II/lac/source/solver_control.cc @@ -81,6 +81,7 @@ SolverControl::check (const unsigned int step, if (_log_result) deallog << "Failure step " << step << " value " << check_value << std::endl; + lcheck = failure; return failure; } @@ -89,14 +90,22 @@ SolverControl::check (const unsigned int step, if (_log_result) deallog << "Convergence step " << step << " value " << check_value << std::endl; + lcheck = success; return success; } - + lcheck = iterate; return iterate; } +SolverControl::State +SolverControl::last_check() const +{ + return lcheck; +} + + double SolverControl::last_value() const { @@ -181,6 +190,7 @@ ReductionControl::check (const unsigned int step, if (_log_result) deallog << "Convergence step " << step << " value " << check_value << std::endl; + lcheck = success; return success; } else