From b8c883a0a74c42888f7bd9c2a8ef39d846c07fec Mon Sep 17 00:00:00 2001
From: Timo Heister <timo.heister@gmail.com>
Date: Tue, 9 Oct 2012 22:08:00 +0000
Subject: [PATCH] treat lucky breakdown in GMRES and FGMRES

git-svn-id: https://svn.dealii.org/trunk@27040 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/doc/news/changes.h                 |  4 ++++
 deal.II/include/deal.II/lac/solver_gmres.h | 11 ++++++++---
 2 files changed, 12 insertions(+), 3 deletions(-)

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 414c76aae9..34da06468c 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -87,6 +87,10 @@ DoFHandler, in particular removal of specializations.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: handle lucky breakdowns in GMRES/FGMRES.
+<br>
+(Bärbel Jansse, Timo Heister, 2012/10/09)
+
 <li> Fixed: GridTools::find_cells_adjacent_to_vertex got into
 trouble with anisotropically refined meshes. This is now fixed.
 <br>
diff --git a/deal.II/include/deal.II/lac/solver_gmres.h b/deal.II/include/deal.II/lac/solver_gmres.h
index 2ca953ca83..0d6c63aa62 100644
--- a/deal.II/include/deal.II/lac/solver_gmres.h
+++ b/deal.II/include/deal.II/lac/solver_gmres.h
@@ -706,8 +706,9 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
 
           const double s = vv.l2_norm();
           h(inner_iteration+1) = s;
-          //s=0 is a lucky breakdown. Handle this somehow decently
-          if(s!=0)
+          //s=0 is a lucky breakdown, the solver will reach convergence,
+	  //but we must not divide by zero here.
+          if (numbers::is_finite(1./s))
             vv *= 1./s;
 
                                            /*  Transformation into
@@ -907,7 +908,11 @@ SolverFGMRES<VECTOR>::solve (
 
       for (unsigned int j=0;j<basis_size;++j)
         {
-          v(j,x).equ(1./a, *aux);
+	  if (numbers::is_finite(1./a)) // treat lucky breakdown
+	    v(j,x).equ(1./a, *aux);
+	  else
+	    v(j,x) = 0.;
+	  
 
           precondition.vmult(z(j,x), v[j]);
           A.vmult(*aux, z[j]);
-- 
2.39.5