From 3bc780285664978cb7d39837d122e654fad3c4d1 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 7 Mar 2007 16:08:18 +0000 Subject: [PATCH] Fix a bug in bicgstab reported last week by Roger Young. git-svn-id: https://svn.dealii.org/trunk@14545 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.html | 7 ++ deal.II/lac/include/lac/solver_bicgstab.h | 15 +++- tests/lac/Makefile | 3 +- tests/lac/bicgstab_early.cc | 91 +++++++++++++++++++++++ tests/lac/bicgstab_early/cmp/generic | 2 + tests/lac/solver/cmp/generic | 15 ++-- 6 files changed, 123 insertions(+), 10 deletions(-) create mode 100644 tests/lac/bicgstab_early.cc create mode 100644 tests/lac/bicgstab_early/cmp/generic diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index 7a17e56e14..deb41f8b55 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -735,6 +735,13 @@ inconvenience this causes.

lac

    +
  1. Fixed: The SolverBicgstab class + did not handle hitting on the solution early very + gracefully (it threw an exception). This is now fixed. +
    + (Roger Young 2007/03/07) +

    +
  2. Fixed: The SparseDirectMA27 class allows to run the sparse direct solver as a separate program (the detached_ma27 program) rather than as part of the current diff --git a/deal.II/lac/include/lac/solver_bicgstab.h b/deal.II/lac/include/lac/solver_bicgstab.h index a1b40372f9..64f29a331b 100644 --- a/deal.II/lac/include/lac/solver_bicgstab.h +++ b/deal.II/lac/include/lac/solver_bicgstab.h @@ -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::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; diff --git a/tests/lac/Makefile b/tests/lac/Makefile index 201ea8aaf1..a6912fcbd5 100644 --- a/tests/lac/Makefile +++ b/tests/lac/Makefile @@ -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 index 0000000000..d6cff9d8a2 --- /dev/null +++ b/tests/lac/bicgstab_early.cc @@ -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 +#include +#include +#include +#include "testmatrix.h" +#include +#include +#include +#include +#include +#include + + +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 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 rhs(4); + rhs(0) = rhs(2) = 0; + rhs(1) = rhs(3) = 0.0975; + + + SparseILU ilu (sparsity_pattern); + ilu.decompose (M); + + Vector 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 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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/lac/bicgstab_early/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/lac/solver/cmp/generic b/tests/lac/solver/cmp/generic index c234a9c709..47e80636c0 100644 --- a/tests/lac/solver/cmp/generic +++ b/tests/lac/solver/cmp/generic @@ -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 -- 2.39.5