]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
solvers throw exception on failure
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 29 Mar 2001 23:35:15 +0000 (23:35 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 29 Mar 2001 23:35:15 +0000 (23:35 +0000)
git-svn-id: https://svn.dealii.org/trunk@4342 0785d39b-7218-0410-832d-ea1e28bc413d

12 files changed:
deal.II/doc/news/2001/c-3-1.html
deal.II/lac/include/lac/eigen.h
deal.II/lac/include/lac/solver.h
deal.II/lac/include/lac/solver_bicgstab.h
deal.II/lac/include/lac/solver_cg.h
deal.II/lac/include/lac/solver_control.h
deal.II/lac/include/lac/solver_gmres.h
deal.II/lac/include/lac/solver_minres.h
deal.II/lac/include/lac/solver_qmrs.h
deal.II/lac/include/lac/solver_richardson.h
deal.II/lac/include/lac/solver_selector.h
deal.II/lac/source/solver_control.cc

index 41224c6513c5d238f95406fcb51ec7e479faf2f3..917b44fa2648c04af9776125f06380fbc5ba78d8 100644 (file)
@@ -238,6 +238,14 @@ documentation, etc</a>.
 <h3>lac</h3>
 
 <ol>
+  <li> <p> Changed: Solver functions <code class="member">solve</code>
+       are void now. If the solver has not converged within the
+       maximum number of iterations, it throws an exception of type
+       <code class="class">ExcNoConvergence</code>.
+       <br>
+       (GK 2001/03/29)
+       </p>
+
   <li> <p> 
        New: The functions <code class="class">FullMatrix</code>::<code
        class="member">mmult</code> and <code
index 3eb0ce291df9b1392e5053c422976844087c56f7..c64c643b8e66eb51e435c24c49ec8279ab65f1e2 100644 (file)
@@ -139,7 +139,7 @@ class EigenPower : private Solver<VECTOR>
                                      * l2-norm.
                                      */
     template <class MATRIX>
-    typename Solver<VECTOR>::ReturnState
+    void
     solve (double       &value,
           const MATRIX &A,
           VECTOR       &x);
@@ -218,7 +218,7 @@ class EigenInverse : private Solver<VECTOR>
                                      * l2-norm.
                                      */
     template <class MATRIX>
-    typename Solver<VECTOR>::ReturnState
+    void
     solve (double       &value,
           const MATRIX &A,
           VECTOR       &x);
@@ -301,7 +301,7 @@ EigenPower<VECTOR>::~EigenPower ()
 
 template <class VECTOR>
 template <class MATRIX>
-typename Solver<VECTOR>::ReturnState
+void
 EigenPower<VECTOR>::solve (double       &value,
                           const MATRIX &A,
                           VECTOR       &x)
@@ -362,10 +362,9 @@ EigenPower<VECTOR>::solve (double       &value,
 
   deallog.pop();
                                   // Output
-  if (conv == SolverControl::failure)
-    return exceeded;
-  else
-    return success;
+  AssertThrow(control().last_check() == SolverControl::success,
+             typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+                                                       control().last_value()));
 }
 
 //----------------------------------------------------------------------//
@@ -386,7 +385,7 @@ EigenInverse<VECTOR>::~EigenInverse ()
 
 template <class VECTOR>
 template <class MATRIX>
-typename Solver<VECTOR>::ReturnState
+void
 EigenInverse<VECTOR>::solve (double       &value,
                             const MATRIX &A,
                             VECTOR       &x)
@@ -472,10 +471,9 @@ EigenInverse<VECTOR>::solve (double       &value,
   deallog.pop();
 
                                   // Output
-  if (conv == SolverControl::failure)
-    return exceeded;
-  else
-    return success;
+  AssertThrow(control().last_check() == SolverControl::success,
+             typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+                                                       control().last_value()));
 }
 
 #endif
index f76a65bfff5a5b7fecb94bd423e0ff4ef9c4b82c..daa52fa5426b2f8e4320ed2b7d056edf80275112 100644 (file)
@@ -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 Vector = ::Vector<double> >
 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:
 
                                     /**
index 63a282e657b37aae7563d59d1c88f57bddf4d571..9e6676879c16fcdef21e0fbce13eb2114fc6bf49 100644 (file)
@@ -97,7 +97,7 @@ class SolverBicgstab : public Subscriptor, private Solver<VECTOR>
                                      * Solve primal problem only.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    typename Solver<VECTOR>::ReturnState
+    void
     solve (const MATRIX &A,
           VECTOR       &x,
           const VECTOR &b,
@@ -212,7 +212,7 @@ class SolverBicgstab : public Subscriptor, private Solver<VECTOR>
                                      * The iteration loop itself.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    typename Solver<VECTOR>::ReturnState 
+    bool
     iterate(const MATRIX& A, const PRECONDITIONER& precondition);
   
 };
@@ -280,7 +280,7 @@ SolverBicgstab<VECTOR>::print_vectors(const unsigned int,
 
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
+bool
 SolverBicgstab<VECTOR>::iterate(const MATRIX& A,
                                const PRECONDITIONER& precondition)
 {
@@ -312,7 +312,7 @@ SolverBicgstab<VECTOR>::iterate(const MATRIX& A,
 //TODO:[GK] Find better breakdown criterion
 
       if (fabs(alpha) > 1.e10)
-       return typename Solver<VECTOR>::ReturnState(breakdown);
+       return true;
     
       s.equ(1., r, -alpha, v);
       precondition.vmult(z,s);
@@ -322,18 +322,18 @@ SolverBicgstab<VECTOR>::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<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
+void
 SolverBicgstab<VECTOR>::solve(const MATRIX &A,
                              VECTOR       &x,
                              const VECTOR &b,
@@ -354,7 +354,7 @@ SolverBicgstab<VECTOR>::solve(const MATRIX &A,
 
   step = 0;
 
-  typename Solver<VECTOR>::ReturnState state = breakdown;
+  bool state;
   
   do 
     {
@@ -363,7 +363,7 @@ SolverBicgstab<VECTOR>::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<VECTOR>::solve(const MATRIX &A,
   memory.free(Vv);
   
   deallog.pop();
-  return state;
+  
+  AssertThrow(control().last_check() == SolverControl::success,
+             typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+                                                       control().last_value()));
 }
 
 
index 79c4369aca8993585fc182882eaceb0c7abcc5fc..bf81acefc43972db9fe7ec578932b80bf1fa38dd 100644 (file)
@@ -16,6 +16,7 @@
 
 #include <lac/solver.h>
 #include <lac/solver_control.h>
+#include <base/exceptions.h>
 #include <base/logstream.h>
 #include <base/subscriptor.h>
 #include <cmath>
@@ -87,7 +88,7 @@ class SolverCG : public Subscriptor, private Solver<VECTOR>
                                      * Solver method.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    typename Solver<VECTOR>::ReturnState
+    void
     solve (const MATRIX &A,
           VECTOR       &x,
           const VECTOR &b,
@@ -187,7 +188,7 @@ SolverCG<VECTOR>::print_vectors(const unsigned int,
 
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState 
+void
 SolverCG<VECTOR>::solve (const MATRIX &A,
                         VECTOR       &x,
                         const VECTOR &b,
@@ -231,7 +232,7 @@ SolverCG<VECTOR>::solve (const MATRIX &A,
       memory.free(Vz);
       memory.free(VAp);
       deallog.pop();
-      return success;
+      return;
     };
   
   g.scale(-1.);
@@ -281,10 +282,9 @@ SolverCG<VECTOR>::solve (const MATRIX &A,
                                   // Output
   deallog.pop();
  
-  if (conv == SolverControl::failure)
-    return exceeded;
-  else
-    return success;
+  AssertThrow(control().last_check() == SolverControl::success,
+             typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+                                                       control().last_value()));
 };
 
 
index 648d8fef703c32ce5032e788455af132be8be668..839e11e2ff031e6a957080f4e47934893cd89f4b 100644 (file)
@@ -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.
                                      */
index 44764f23ad86b349484a13fa8b5402456b645bef..cef56c00665e01addb2550e5a8a7a994bc04bd68 100644 (file)
@@ -110,10 +110,10 @@ class SolverGMRES : public Subscriptor, private Solver<VECTOR>
                                      * Solver method.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    typename Solver<VECTOR>::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<VECTOR>::givens_rotation (Vector<double> &h,
 
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
-SolverGMRES<VECTOR>::solve (const MATRIX         &A,
-                           VECTOR               &x,
-                           const VECTOR         &b,
-                           const PRECONDITIONER &precondition)
+void
+SolverGMRES<VECTOR>::solve (const MATRIXA,
+                                  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<VECTOR>::ExcNoConvergence(control().last_step(),
+                                                       control().last_value()));
 };
 
 
index 163121e4b33e21a79b906ae37bc4a97fd0cba2aa..fa7f86742c70bbe17a9a94435d3734fd3caa1cbf 100644 (file)
@@ -77,7 +77,7 @@ class SolverMinRes : public Subscriptor, private Solver<VECTOR>
                                      * Solver method.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    typename Solver<VECTOR>::ReturnState
+    void
     solve (const MATRIX &A,
           VECTOR       &x,
           const VECTOR &b,
@@ -351,10 +351,9 @@ SolverMinRes<VECTOR>::solve (const MATRIX &A,
                                   // Output
   deallog.pop ();
   
-  if (conv == SolverControl::failure)
-    return exceeded;
-  
-  return success;
+  AssertThrow(control().last_check() == SolverControl::success,
+             typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+                                                       control().last_value()));
 };
 
 
index fa6f5d1c6e1703cc4cf349251c21bebce0dc2cc1..4aa2aedf9bb44a258041533e2d16804b7edf063d 100644 (file)
@@ -109,11 +109,11 @@ class SolverQMRS : public Subscriptor, private Solver<VECTOR>
                                      * Solver method.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    typename Solver<VECTOR>::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<VECTOR>
                                      * The iteration loop itself.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    typename Solver<VECTOR>::ReturnState 
+    bool
     iterate(const MATRIX& A, const PRECONDITIONER& precondition);
                                     /**
                                      * The current iteration step.
@@ -217,7 +217,7 @@ SolverQMRS<VECTOR>::print_vectors(const unsigned int,
 
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState 
+void
 SolverQMRS<VECTOR>::solve (const MATRIX &A,
                           VECTOR       &x,
                           const VECTOR &b,
@@ -244,7 +244,7 @@ SolverQMRS<VECTOR>::solve (const MATRIX &A,
 
   step = 0;
   
-  typename Solver<VECTOR>::ReturnState state = breakdown;
+  bool state;
 
   do 
     {
@@ -252,7 +252,7 @@ SolverQMRS<VECTOR>::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<VECTOR>::solve (const MATRIX &A,
   // Output
 
   deallog.pop();
-  return state;
+
+  AssertThrow(control().last_check() == SolverControl::success,
+             typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+                                                       control().last_value()));
 };
 
 
 
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState 
+bool
 SolverQMRS<VECTOR>::iterate(const MATRIX& A,
                            const PRECONDITIONER& precondition)
 {
@@ -309,7 +311,7 @@ SolverQMRS<VECTOR>::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;
 }
 
 
index 847f2b40971a88eda4cab24fa55ade0f4e701ee5..103550ec44dc91ee86ab3df0066e576623b6d6a8 100644 (file)
@@ -79,19 +79,19 @@ class SolverRichardson : public Subscriptor, private Solver<VECTOR>
                                      * Solve $Ax=b$ for $x$.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    typename Solver<VECTOR>::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<class MATRIX, class PRECONDITIONER>
-    typename Solver<VECTOR>::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<VECTOR>::~SolverRichardson()
 
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState 
+void
 SolverRichardson<VECTOR>::solve (const MATRIX &A,
                                 VECTOR       &x,
                                 const VECTOR &b,
@@ -206,16 +206,15 @@ SolverRichardson<VECTOR>::solve (const MATRIX &A,
 
   deallog.pop();
                                   // Output
-  if (conv == SolverControl::failure)
-    return exceeded;
-  else
-    return success;
+  AssertThrow(control().last_check() == SolverControl::success,
+             typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+                                                       control().last_value()));
 }
 
 
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState 
+void
 SolverRichardson<VECTOR>::Tsolve (const MATRIX &A,
                                  VECTOR       &x,
                                  const VECTOR &b,
@@ -253,10 +252,9 @@ SolverRichardson<VECTOR>::Tsolve (const MATRIX &A,
 
   deallog.pop();
                                   // Output
-  if (conv == SolverControl::failure)
-    return exceeded;
-  else
-    return success;
+  AssertThrow(control().last_check() == SolverControl::success,
+             typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+                                                       control().last_value()));
 }
 
 
index be805096f7360e009f75bdf075d1912b2807e5a2..268c441e942dccf758e83b02f4e2a328e970719c 100644 (file)
@@ -111,10 +111,10 @@ class SolverSelector : public Subscriptor
                                      * 
                                      */
     template<class Matrix, class Preconditioner>
-    typename Solver<Vector>::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<Vector>::~SolverSelector()
 
 template <class Vector>
 template<class Matrix, class Preconditioner>
-typename Solver<Vector>::ReturnState 
+void
 SolverSelector<Vector>::solve(const Matrix &A,
                              Vector &x,
                              const Vector &b,
@@ -230,30 +230,28 @@ SolverSelector<Vector>::solve(const Matrix &A,
     {
       SolverRichardson<Vector> 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<Vector> 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<Vector> 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<Vector> 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<Vector>::breakdown;
 };
 
 
index e9c2e9c6c48cff053e06ee9821a571ac409b3edf..27f671efb7e02ff84a0dba3d7102a1b3d84337f6 100644 (file)
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.