]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix a bug in bicgstab reported last week by Roger Young.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 7 Mar 2007 16:08:18 +0000 (16:08 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 7 Mar 2007 16:08:18 +0000 (16:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@14545 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.html
deal.II/lac/include/lac/solver_bicgstab.h
tests/lac/Makefile
tests/lac/bicgstab_early.cc [new file with mode: 0644]
tests/lac/bicgstab_early/cmp/generic [new file with mode: 0644]
tests/lac/solver/cmp/generic

index 7a17e56e14fe02f3df17d7b5177a5b2d5263c517..deb41f8b5570eab31624ef8ea6cca7348b4d7e1f 100644 (file)
@@ -735,6 +735,13 @@ inconvenience this causes.
 <h3>lac</h3>
 
 <ol>
+  <li> <p>Fixed: The <code class="class">SolverBicgstab</code> class
+       did not handle hitting on the solution early very
+       gracefully (it threw an exception). This is now fixed.
+       <br>
+       (Roger Young 2007/03/07)
+       </p>
+
   <li> <p>Fixed: The <code class="member">SparseDirectMA27</code> class allows
        to run the sparse direct solver as a separate program (the
        <code>detached_ma27</code> program) rather than as part of the current
index a1b40372f988468e72b482ad4d0e17f2d61edb5e..64f29a331ba12f63ac26e39e9fe04d389281bfe5 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -352,6 +352,19 @@ SolverBicgstab<VECTOR>::iterate(const MATRIX& A,
        return true;
       
       s.equ(1., r, -alpha, v);
+
+                                      // check for early success, see
+                                      // the lac/bicgstab_early
+                                      // testcase as to why this is
+                                      // necessary
+      if (this->control().check(step, s.l2_norm()/Vb->l2_norm())
+         == SolverControl::success)
+       {
+         Vx->add(alpha, y);
+         print_vectors(step, *Vx, r, y);
+         return false;
+       }
+
       precondition.vmult(z,s);
       A.vmult(t,z);
       rhobar = t*s;
index 201ea8aaf1b69202b61d47a7c733f6d524c79008..a6912fcbd51ffaa37d15bd9075e9f63622659d1f 100644 (file)
@@ -30,7 +30,8 @@ tests_x = vector-vector \
        solver eigen \
        sparse_ilu sparse_ilu_t lapack_fill block_minres \
        identity_matrix* trace \
-       print_formatted_*
+       print_formatted_* \
+       bicgstab_early
 
 
 # from above list of regular expressions, generate the real set of
diff --git a/tests/lac/bicgstab_early.cc b/tests/lac/bicgstab_early.cc
new file mode 100644 (file)
index 0000000..d6cff9d
--- /dev/null
@@ -0,0 +1,91 @@
+//----------------------------  bicgstab_early.cc  ---------------------------
+//    $Id: bicgstab_early.cc 11749 2005-11-09 19:11:20Z wolf $
+//    Version: $Name$ 
+//
+//    Copyright (C) 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  bicgstab_early.cc  ---------------------------
+
+// adapted from a testcase by Roger Young, sent to the mailing list
+// 2007-03-02, that illustrates that bicgstab can't handle early
+// success
+
+#include "../tests.h"
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include "testmatrix.h"
+#include <base/logstream.h>
+#include <lac/sparse_matrix.h>
+#include <lac/sparse_ilu.h>
+#include <lac/solver_bicgstab.h>
+#include <lac/vector.h>
+#include <lac/precondition.h>
+
+
+int main()
+{
+  std::ofstream logfile("bicgstab_early/output");
+  logfile.precision(4);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+  
+  GrowingVectorMemory<> mem;
+  SolverControl control(100, 1.e-3);
+  SolverBicgstab<> bicgstab(control, mem);
+
+  SparsityPattern sparsity_pattern(4,4,4);
+  for (unsigned int i=0; i<4; ++i)
+    for (unsigned int j=0; j<4; ++j)
+      sparsity_pattern.add(i,j);
+  sparsity_pattern.compress();
+  
+  SparseMatrix<double> M(sparsity_pattern);
+  M.set(0,0,21.1);
+  M.set(0,1,0);
+  M.set(0,2,0);
+  M.set(0,3,0);
+  M.set(1,1,7.033333333);
+  M.set(1,0,0);
+  M.set(1,2,0);
+  M.set(1,3,3.516666667);
+  M.set(2,2,21.1);
+  M.set(2,0,0);
+  M.set(2,1,0);
+  M.set(2,3,0);
+  M.set(3,3,7.033333333);
+  M.set(3,0,0);
+  M.set(3,1,3.516666667);
+  M.set(3,2,0);
+
+  Vector<double> rhs(4);
+  rhs(0) = rhs(2) = 0;
+  rhs(1) = rhs(3) = 0.0975;
+
+
+  SparseILU<double> ilu (sparsity_pattern);
+  ilu.decompose (M);
+
+  Vector<double> solution (4);
+                                  // this would fail with elements of
+                                  // the solution vector being set to
+                                  // NaN before Roger's suggested
+                                  // change
+  bicgstab.solve (M, solution, rhs, ilu);
+
+  for (unsigned int i=0; i<4; ++i)
+    deallog << solution(i) << std::endl;
+
+  Vector<double> res (4);
+  M.residual (res, solution, rhs);
+  deallog << "residual=" << res.l2_norm()
+         << std::endl;
+}
+
diff --git a/tests/lac/bicgstab_early/cmp/generic b/tests/lac/bicgstab_early/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
index c234a9c70964cdf6e80d5f7dc6086acbfd36e0d6..47e80636c0630e2542708edc62d79b5d3c5b0390 100644 (file)
@@ -70,7 +70,7 @@ DEAL:psor:cg::Starting value 3.000
 DEAL:psor:cg::Failure step 100 value 0.1024
 DEAL:psor::Iterative method reported convergence failure in step 100 with residual 0.102391
 DEAL:psor:Bicgstab::Starting value 3.000
-DEAL:psor:Bicgstab::Convergence step 4 value 0.0001279
+DEAL:psor:Bicgstab::Convergence step 4 value 0.0002656
 DEAL:psor:GMRES::Starting value 1.467
 DEAL:psor:GMRES::Convergence step 5 value 0.0006649
 DEAL:psor:GMRES::Starting value 1.467
@@ -81,8 +81,7 @@ DEAL:no-fail:cg::Starting value 11.00
 DEAL:no-fail:cg::Failure step 10 value 0.1496
 DEAL:no-fail::Iterative method reported convergence failure in step 10 with residual 0.149566
 DEAL:no-fail:Bicgstab::Starting value 11.00
-DEAL:no-fail:Bicgstab::Failure step 10 value 0.001961
-DEAL:no-fail::Iterative method reported convergence failure in step 10 with residual 0.00196051
+DEAL:no-fail:Bicgstab::Convergence step 10 value 0.0002572
 DEAL:no-fail:GMRES::Starting value 11.00
 DEAL:no-fail:GMRES::Failure step 10 value 0.8414
 DEAL:no-fail::Iterative method reported convergence failure in step 10 with residual 0.841354
@@ -98,7 +97,7 @@ DEAL:no::Iterative method reported convergence failure in step 100 with residual
 DEAL:no:cg::Starting value 11.00
 DEAL:no:cg::Convergence step 15 value 0.0005794
 DEAL:no:Bicgstab::Starting value 11.00
-DEAL:no:Bicgstab::Convergence step 11 value 0.0003990
+DEAL:no:Bicgstab::Convergence step 10 value 0.0002572
 DEAL:no:GMRES::Starting value 11.00
 DEAL:no:GMRES::Convergence step 43 value 0.0009788
 DEAL:no:GMRES::Starting value 11.00
@@ -111,7 +110,7 @@ DEAL:rich::Iterative method reported convergence failure in step 100 with residu
 DEAL:rich:cg::Starting value 11.00
 DEAL:rich:cg::Convergence step 15 value 0.0005794
 DEAL:rich:Bicgstab::Starting value 11.00
-DEAL:rich:Bicgstab::Convergence step 11 value 0.0003990
+DEAL:rich:Bicgstab::Convergence step 10 value 0.0002572
 DEAL:rich:GMRES::Starting value 6.600
 DEAL:rich:GMRES::Convergence step 42 value 0.0007803
 DEAL:rich:GMRES::Starting value 6.600
@@ -125,7 +124,7 @@ DEAL:ssor:Richardson::Convergence step 59 value 0.0008892
 DEAL:ssor:cg::Starting value 11.00
 DEAL:ssor:cg::Convergence step 8 value 0.0005816
 DEAL:ssor:Bicgstab::Starting value 11.00
-DEAL:ssor:Bicgstab::Convergence step 5 value 9.315e-05
+DEAL:ssor:Bicgstab::Convergence step 4 value 0.0003834
 DEAL:ssor:GMRES::Starting value 11.05
 DEAL:ssor:GMRES::Convergence step 7 value 0.0006838
 DEAL:ssor:GMRES::Starting value 11.05
@@ -140,7 +139,7 @@ DEAL:sor:cg::Starting value 11.00
 DEAL:sor:cg::Failure step 100 value 5.643
 DEAL:sor::Iterative method reported convergence failure in step 100 with residual 5.64329
 DEAL:sor:Bicgstab::Starting value 11.00
-DEAL:sor:Bicgstab::Convergence step 14 value 0.0009987
+DEAL:sor:Bicgstab::Convergence step 12 value 0.0007246
 DEAL:sor:GMRES::Starting value 7.322
 DEAL:sor:GMRES::Convergence step 23 value 0.0006981
 DEAL:sor:GMRES::Starting value 7.322
@@ -153,7 +152,7 @@ DEAL:psor:cg::Starting value 11.00
 DEAL:psor:cg::Failure step 100 value 2.935
 DEAL:psor::Iterative method reported convergence failure in step 100 with residual 2.93488
 DEAL:psor:Bicgstab::Starting value 11.00
-DEAL:psor:Bicgstab::Convergence step 11 value 0.0001588
+DEAL:psor:Bicgstab::Convergence step 10 value 0.0003167
 DEAL:psor:GMRES::Starting value 7.345
 DEAL:psor:GMRES::Convergence step 20 value 0.0008491
 DEAL:psor:GMRES::Starting value 7.345

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.