From: Peter Munch Date: Sat, 12 Nov 2022 15:37:44 +0000 (+0100) Subject: Fix vector access in SolverGMRES X-Git-Tag: v9.5.0-rc1~727^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a0a8446705a9ee336e46f7915bf9a92a5012db5d;p=dealii.git Fix vector access in SolverGMRES --- diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h index 76bc111c58..e80704e173 100644 --- a/include/deal.II/lac/solver_gmres.h +++ b/include/deal.II/lac/solver_gmres.h @@ -818,10 +818,10 @@ namespace internal { Assert(dim > 0, ExcInternalError()); - for (unsigned int i = 0; i < dim - 1; ++i) + for (unsigned int i = 0; i < dim; ++i) vv.add(-h(i), orthogonal_vectors[i]); - return vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv); + return std::sqrt(vv.add_and_dot(-h(dim), orthogonal_vectors[dim], vv)); } @@ -891,10 +891,10 @@ namespace internal for (; j < vv.locally_owned_size(); ++j) { - double temp = vv(j); + double temp = vv.local_element(j); for (unsigned int i = 0; i < dim; ++i) - temp -= h(i) * orthogonal_vectors[i](j); - vv(j) = temp; + temp -= h(i) * orthogonal_vectors[i].local_element(j); + vv.local_element(j) = temp; norm_vv_temp += temp * temp; } diff --git a/tests/lac/solver.cc b/tests/lac/solver.cc index 354fb8678a..0ddfc6c3d3 100644 --- a/tests/lac/solver.cc +++ b/tests/lac/solver.cc @@ -100,6 +100,11 @@ main() SolverQMRS<> qmrs(control, mem); SolverFIRE<> fire(control, mem); + SolverGMRES<>::AdditionalData data3(8); + data3.orthogonalization_strategy = SolverGMRES<>::AdditionalData:: + OrthogonalizationStrategy::classical_gram_schmidt; + SolverGMRES<> gmresclassical(control, mem, data3); + for (unsigned int size = 4; size <= 30; size *= 3) { unsigned int dim = (size - 1) * (size - 1); @@ -169,6 +174,7 @@ main() check_solve(bicgstab, A, u, f, prec_no); check_solve(gmres, A, u, f, prec_no); check_solve(gmresright, A, u, f, prec_no); + check_solve(gmresclassical, A, u, f, prec_no); // check_solve(minres,A,u,f,prec_no); check_solve(qmrs, A, u, f, prec_no); @@ -187,6 +193,7 @@ main() check_solve(bicgstab, A, u, f, prec_no); check_solve(gmres, A, u, f, prec_no); check_solve(gmresright, A, u, f, prec_no); + check_solve(gmresclassical, A, u, f, prec_no); check_solve(qmrs, A, u, f, prec_no); check_solve(fire, A, u, f, prec_no); rich.set_omega(1.); @@ -201,6 +208,7 @@ main() check_solve(bicgstab, A, u, f, prec_richardson); check_solve(gmres, A, u, f, prec_richardson); check_solve(gmresright, A, u, f, prec_richardson); + check_solve(gmresclassical, A, u, f, prec_richardson); check_solve(qmrs, A, u, f, prec_richardson); check_solve(fire, A, u, f, prec_richardson); rich.set_omega(1.); @@ -215,6 +223,7 @@ main() check_solve(bicgstab, A, u, f, prec_ssor); check_solve(gmres, A, u, f, prec_ssor); check_solve(gmresright, A, u, f, prec_ssor); + check_solve(gmresclassical, A, u, f, prec_ssor); check_solve(qmrs, A, u, f, prec_ssor); check_solve(fire, A, u, f, prec_ssor); @@ -228,6 +237,7 @@ main() check_solve(bicgstab, A, u, f, prec_sor); check_solve(gmres, A, u, f, prec_sor); check_solve(gmresright, A, u, f, prec_sor); + check_solve(gmresclassical, A, u, f, prec_sor); check_solve(fire, A, u, f, prec_sor); deallog.pop(); @@ -240,6 +250,7 @@ main() check_solve(bicgstab, A, u, f, prec_psor); check_solve(gmres, A, u, f, prec_psor); check_solve(gmresright, A, u, f, prec_psor); + check_solve(gmresclassical, A, u, f, prec_psor); check_solve(fire, A, u, f, prec_psor); deallog.pop(); diff --git a/tests/lac/solver.output b/tests/lac/solver.output index cacd06ee37..653812ce53 100644 --- a/tests/lac/solver.output +++ b/tests/lac/solver.output @@ -2,43 +2,49 @@ DEAL::Size 4 Unknowns 9 DEAL::SOR-diff:0.000 DEAL:no-fail:cg::Starting value 3.000 -DEAL:no-fail:cg::Convergence step 3 value 5.917e-17 +DEAL:no-fail:cg::Convergence step 3 value 4.807e-17 DEAL:no-fail:Bicgstab::Starting value 3.000 -DEAL:no-fail:Bicgstab::Convergence step 3 value 4.768e-18 +DEAL:no-fail:Bicgstab::Convergence step 3 value 1.104e-17 DEAL:no-fail:GMRES::Starting value 3.000 -DEAL:no-fail:GMRES::Convergence step 3 value 6.160e-16 +DEAL:no-fail:GMRES::Convergence step 3 value 4.551e-16 DEAL:no-fail:GMRES::Starting value 3.000 -DEAL:no-fail:GMRES::Convergence step 3 value 6.160e-16 +DEAL:no-fail:GMRES::Convergence step 3 value 4.551e-16 +DEAL:no-fail:GMRES::Starting value 3.000 +DEAL:no-fail:GMRES::Convergence step 3 value 4.414e-16 DEAL:no-fail:SQMR::Starting value 3.000 -DEAL:no-fail:SQMR::Convergence step 3 value 9.088e-16 +DEAL:no-fail:SQMR::Convergence step 3 value 1.280e-15 DEAL:no-fail:FIRE::Starting value 9.000 DEAL:no-fail:FIRE::Convergence step 38 value 0.0002184 DEAL:no:Richardson::Starting value 3.000 DEAL:no:Richardson::Convergence step 24 value 0.0007118 DEAL:no:cg::Starting value 3.000 -DEAL:no:cg::Convergence step 3 value 5.917e-17 +DEAL:no:cg::Convergence step 3 value 4.807e-17 DEAL:no:Bicgstab::Starting value 3.000 -DEAL:no:Bicgstab::Convergence step 3 value 4.768e-18 +DEAL:no:Bicgstab::Convergence step 3 value 1.104e-17 +DEAL:no:GMRES::Starting value 3.000 +DEAL:no:GMRES::Convergence step 3 value 4.551e-16 DEAL:no:GMRES::Starting value 3.000 -DEAL:no:GMRES::Convergence step 3 value 6.160e-16 +DEAL:no:GMRES::Convergence step 3 value 4.551e-16 DEAL:no:GMRES::Starting value 3.000 -DEAL:no:GMRES::Convergence step 3 value 6.160e-16 +DEAL:no:GMRES::Convergence step 3 value 4.414e-16 DEAL:no:SQMR::Starting value 3.000 -DEAL:no:SQMR::Convergence step 3 value 9.088e-16 +DEAL:no:SQMR::Convergence step 3 value 1.280e-15 DEAL:no:FIRE::Starting value 9.000 DEAL:no:FIRE::Convergence step 38 value 0.0002184 DEAL:rich:Richardson::Starting value 3.000 DEAL:rich:Richardson::Convergence step 42 value 0.0008696 DEAL:rich:cg::Starting value 3.000 -DEAL:rich:cg::Convergence step 3 value 2.009e-16 +DEAL:rich:cg::Convergence step 3 value 2.871e-16 DEAL:rich:Bicgstab::Starting value 3.000 -DEAL:rich:Bicgstab::Convergence step 3 value 1.708e-17 +DEAL:rich:Bicgstab::Convergence step 3 value 3.238e-18 DEAL:rich:GMRES::Starting value 1.800 -DEAL:rich:GMRES::Convergence step 3 value 2.927e-16 +DEAL:rich:GMRES::Convergence step 3 value 3.294e-16 DEAL:rich:GMRES::Starting value 3.000 -DEAL:rich:GMRES::Convergence step 3 value 6.449e-16 +DEAL:rich:GMRES::Convergence step 3 value 7.148e-16 +DEAL:rich:GMRES::Starting value 1.800 +DEAL:rich:GMRES::Convergence step 3 value 1.934e-16 DEAL:rich:SQMR::Starting value 3.000 -DEAL:rich:SQMR::Convergence step 3 value 1.351e-15 +DEAL:rich:SQMR::Convergence step 3 value 1.422e-15 DEAL:rich:FIRE::Starting value 9.000 DEAL:rich:FIRE::Convergence step 32 value 0.0009708 DEAL:ssor:RichardsonT::Starting value 3.000 @@ -53,6 +59,8 @@ DEAL:ssor:GMRES::Starting value 1.920 DEAL:ssor:GMRES::Convergence step 3 value 0.0003574 DEAL:ssor:GMRES::Starting value 3.000 DEAL:ssor:GMRES::Convergence step 4 value 5.994e-05 +DEAL:ssor:GMRES::Starting value 1.920 +DEAL:ssor:GMRES::Convergence step 3 value 0.0003574 DEAL:ssor:SQMR::Starting value 3.000 DEAL:ssor:SQMR::Convergence step 4 value 0.0001345 DEAL:ssor:FIRE::Starting value 9.000 @@ -65,11 +73,13 @@ DEAL:sor:cg::Starting value 3.000 DEAL:sor:cg::Failure step 100 value 0.2585 DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:sor:Bicgstab::Starting value 3.000 -DEAL:sor:Bicgstab::Convergence step 5 value 4.201e-18 +DEAL:sor:Bicgstab::Convergence step 5 value 4.331e-18 DEAL:sor:GMRES::Starting value 1.462 -DEAL:sor:GMRES::Convergence step 5 value 2.605e-18 +DEAL:sor:GMRES::Convergence step 5 value 9.441e-19 DEAL:sor:GMRES::Starting value 3.000 -DEAL:sor:GMRES::Convergence step 5 value 5.190e-18 +DEAL:sor:GMRES::Convergence step 5 value 1.393e-17 +DEAL:sor:GMRES::Starting value 1.462 +DEAL:sor:GMRES::Convergence step 21 value 0.0007019 DEAL:sor:FIRE::Starting value 9.000 DEAL:sor:FIRE::Convergence step 51 value 0.0004993 DEAL:psor:RichardsonT::Starting value 3.000 @@ -85,6 +95,8 @@ DEAL:psor:GMRES::Starting value 1.467 DEAL:psor:GMRES::Convergence step 5 value 0.0006649 DEAL:psor:GMRES::Starting value 3.000 DEAL:psor:GMRES::Convergence step 6 value 0.0004455 +DEAL:psor:GMRES::Starting value 1.467 +DEAL:psor:GMRES::Convergence step 5 value 0.0006649 DEAL:psor:FIRE::Starting value 9.000 DEAL:psor:FIRE::Convergence step 45 value 0.0009217 DEAL::Size 12 Unknowns 121 @@ -102,6 +114,9 @@ DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, la DEAL:no-fail:GMRES::Starting value 11.00 DEAL:no-fail:GMRES::Failure step 10 value 0.8414 DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res) +DEAL:no-fail:GMRES::Starting value 11.00 +DEAL:no-fail:GMRES::Failure step 10 value 0.8414 +DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res) DEAL:no-fail:SQMR::Starting value 11.00 DEAL:no-fail:SQMR::Failure step 10 value 0.4215 DEAL:no-fail::Exception: SolverControl::NoConvergence(step, state.last_residual) @@ -119,6 +134,8 @@ DEAL:no:GMRES::Starting value 11.00 DEAL:no:GMRES::Convergence step 43 value 0.0009788 DEAL:no:GMRES::Starting value 11.00 DEAL:no:GMRES::Convergence step 43 value 0.0009788 +DEAL:no:GMRES::Starting value 11.00 +DEAL:no:GMRES::Convergence step 43 value 0.0009788 DEAL:no:SQMR::Starting value 11.00 DEAL:no:SQMR::Convergence step 16 value 0.0002583 DEAL:no:FIRE::Starting value 121.0 @@ -135,6 +152,8 @@ DEAL:rich:GMRES::Starting value 6.600 DEAL:rich:GMRES::Convergence step 42 value 0.0007803 DEAL:rich:GMRES::Starting value 11.00 DEAL:rich:GMRES::Convergence step 43 value 0.0009788 +DEAL:rich:GMRES::Starting value 6.600 +DEAL:rich:GMRES::Convergence step 42 value 0.0007803 DEAL:rich:SQMR::Starting value 11.00 DEAL:rich:SQMR::Convergence step 16 value 0.0002583 DEAL:rich:FIRE::Starting value 121.0 @@ -151,6 +170,8 @@ DEAL:ssor:GMRES::Starting value 13.27 DEAL:ssor:GMRES::Convergence step 7 value 0.0008206 DEAL:ssor:GMRES::Starting value 11.00 DEAL:ssor:GMRES::Convergence step 8 value 0.0004929 +DEAL:ssor:GMRES::Starting value 13.27 +DEAL:ssor:GMRES::Convergence step 7 value 0.0008206 DEAL:ssor:SQMR::Starting value 11.00 DEAL:ssor:SQMR::Convergence step 9 value 0.0004140 DEAL:ssor:FIRE::Starting value 121.0 @@ -168,6 +189,8 @@ DEAL:sor:GMRES::Starting value 7.322 DEAL:sor:GMRES::Convergence step 23 value 0.0006981 DEAL:sor:GMRES::Starting value 11.00 DEAL:sor:GMRES::Convergence step 24 value 0.0007120 +DEAL:sor:GMRES::Starting value 7.322 +DEAL:sor:GMRES::Convergence step 23 value 0.0006981 DEAL:sor:FIRE::Starting value 121.0 DEAL:sor:FIRE::Convergence step 92 value 0.0009080 DEAL:psor:RichardsonT::Starting value 11.00 @@ -183,5 +206,7 @@ DEAL:psor:GMRES::Starting value 7.345 DEAL:psor:GMRES::Convergence step 20 value 0.0008491 DEAL:psor:GMRES::Starting value 11.00 DEAL:psor:GMRES::Convergence step 23 value 0.0007600 +DEAL:psor:GMRES::Starting value 7.345 +DEAL:psor:GMRES::Convergence step 20 value 0.0008491 DEAL:psor:FIRE::Starting value 121.0 DEAL:psor:FIRE::Convergence step 93 value 0.0006976 diff --git a/tests/lac/solver_gmres_01.cc b/tests/lac/solver_gmres_01.cc new file mode 100644 index 0000000000..586982087d --- /dev/null +++ b/tests/lac/solver_gmres_01.cc @@ -0,0 +1,92 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Check the path of SolverGMRES with distributed vectors and +// OrthogonalizationStrategy::classical_gram_schmidt. + + +#include +#include +#include +#include + +#include "../tests.h" + + +struct MyDiagonalMatrix +{ + MyDiagonalMatrix(const LinearAlgebra::distributed::Vector &diagonal) + : diagonal(diagonal) + {} + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + dst = src; + dst.scale(diagonal); + } + + const LinearAlgebra::distributed::Vector &diagonal; +}; + + + +SolverControl::State +monitor_norm(const unsigned int iteration, + const double check_value, + const LinearAlgebra::distributed::Vector &) +{ + deallog << " estimated residual at iteration " << iteration << ": " + << check_value << std::endl; + return SolverControl::success; +} + + +int +main() +{ + initlog(); + + // Create diagonal matrix with entries between 1 and 30 + DiagonalMatrix> unit_matrix; + unit_matrix.get_vector().reinit(30); + unit_matrix.get_vector() = 1.0; + + LinearAlgebra::distributed::Vector matrix_entries(unit_matrix.m()); + for (unsigned int i = 0; i < unit_matrix.m(); ++i) + matrix_entries(i) = i + 1; + MyDiagonalMatrix matrix(matrix_entries); + + LinearAlgebra::distributed::Vector rhs(unit_matrix.m()), + sol(unit_matrix.m()); + rhs = 1.; + + deallog << "Solve with PreconditionIdentity: " << std::endl; + SolverControl control(40, 1e-4); + SolverGMRES>::AdditionalData data3( + 8); + data3.orthogonalization_strategy = + SolverGMRES>::AdditionalData:: + OrthogonalizationStrategy::classical_gram_schmidt; + SolverGMRES> solver(control, + data3); + solver.connect(&monitor_norm); + solver.solve(matrix, sol, rhs, PreconditionIdentity()); + + deallog << "Solve with diagonal preconditioner: " << std::endl; + sol = 0; + solver.solve(matrix, sol, rhs, unit_matrix); +} diff --git a/tests/lac/solver_gmres_01.output b/tests/lac/solver_gmres_01.output new file mode 100644 index 0000000000..fcf6589793 --- /dev/null +++ b/tests/lac/solver_gmres_01.output @@ -0,0 +1,89 @@ + +DEAL::Solve with PreconditionIdentity: +DEAL:GMRES::Starting value 5.47723 +DEAL:GMRES:: estimated residual at iteration 0: 5.47723 +DEAL:GMRES:: estimated residual at iteration 1: 2.67042 +DEAL:GMRES:: estimated residual at iteration 2: 1.70538 +DEAL:GMRES:: estimated residual at iteration 3: 1.20190 +DEAL:GMRES:: estimated residual at iteration 4: 0.884567 +DEAL:GMRES:: estimated residual at iteration 5: 0.662215 +DEAL:GMRES:: estimated residual at iteration 6: 0.496462 +DEAL:GMRES:: estimated residual at iteration 6: 0.496462 +DEAL:GMRES:: estimated residual at iteration 7: 0.419909 +DEAL:GMRES:: estimated residual at iteration 8: 0.339229 +DEAL:GMRES:: estimated residual at iteration 9: 0.260767 +DEAL:GMRES:: estimated residual at iteration 10: 0.189115 +DEAL:GMRES:: estimated residual at iteration 11: 0.126467 +DEAL:GMRES:: estimated residual at iteration 12: 0.0736025 +DEAL:GMRES:: estimated residual at iteration 12: 0.0736025 +DEAL:GMRES:: estimated residual at iteration 13: 0.0505275 +DEAL:GMRES:: estimated residual at iteration 14: 0.0360998 +DEAL:GMRES:: estimated residual at iteration 15: 0.0253697 +DEAL:GMRES:: estimated residual at iteration 16: 0.0184200 +DEAL:GMRES:: estimated residual at iteration 17: 0.0142333 +DEAL:GMRES:: estimated residual at iteration 18: 0.0116662 +DEAL:GMRES:: estimated residual at iteration 18: 0.0116662 +DEAL:GMRES:: estimated residual at iteration 19: 0.00989870 +DEAL:GMRES:: estimated residual at iteration 20: 0.00800703 +DEAL:GMRES:: estimated residual at iteration 21: 0.00604107 +DEAL:GMRES:: estimated residual at iteration 22: 0.00431736 +DEAL:GMRES:: estimated residual at iteration 23: 0.00304828 +DEAL:GMRES:: estimated residual at iteration 24: 0.00203679 +DEAL:GMRES:: estimated residual at iteration 24: 0.00203679 +DEAL:GMRES:: estimated residual at iteration 25: 0.00141669 +DEAL:GMRES:: estimated residual at iteration 26: 0.00101755 +DEAL:GMRES:: estimated residual at iteration 27: 0.000724887 +DEAL:GMRES:: estimated residual at iteration 28: 0.000535249 +DEAL:GMRES:: estimated residual at iteration 29: 0.000423665 +DEAL:GMRES:: estimated residual at iteration 30: 0.000358414 +DEAL:GMRES:: estimated residual at iteration 30: 0.000358414 +DEAL:GMRES:: estimated residual at iteration 31: 0.000307521 +DEAL:GMRES:: estimated residual at iteration 32: 0.000247031 +DEAL:GMRES:: estimated residual at iteration 33: 0.000183614 +DEAL:GMRES:: estimated residual at iteration 34: 0.000129847 +DEAL:GMRES::Convergence step 35 value 9.20930e-05 +DEAL:GMRES:: estimated residual at iteration 35: 9.20930e-05 +DEAL::Solve with diagonal preconditioner: +DEAL:GMRES::Starting value 5.47723 +DEAL:GMRES:: estimated residual at iteration 0: 5.47723 +DEAL:GMRES:: estimated residual at iteration 1: 2.67042 +DEAL:GMRES:: estimated residual at iteration 2: 1.70538 +DEAL:GMRES:: estimated residual at iteration 3: 1.20190 +DEAL:GMRES:: estimated residual at iteration 4: 0.884567 +DEAL:GMRES:: estimated residual at iteration 5: 0.662215 +DEAL:GMRES:: estimated residual at iteration 6: 0.496462 +DEAL:GMRES:: estimated residual at iteration 6: 0.496462 +DEAL:GMRES:: estimated residual at iteration 7: 0.419909 +DEAL:GMRES:: estimated residual at iteration 8: 0.339229 +DEAL:GMRES:: estimated residual at iteration 9: 0.260767 +DEAL:GMRES:: estimated residual at iteration 10: 0.189115 +DEAL:GMRES:: estimated residual at iteration 11: 0.126467 +DEAL:GMRES:: estimated residual at iteration 12: 0.0736025 +DEAL:GMRES:: estimated residual at iteration 12: 0.0736025 +DEAL:GMRES:: estimated residual at iteration 13: 0.0505275 +DEAL:GMRES:: estimated residual at iteration 14: 0.0360998 +DEAL:GMRES:: estimated residual at iteration 15: 0.0253697 +DEAL:GMRES:: estimated residual at iteration 16: 0.0184200 +DEAL:GMRES:: estimated residual at iteration 17: 0.0142333 +DEAL:GMRES:: estimated residual at iteration 18: 0.0116662 +DEAL:GMRES:: estimated residual at iteration 18: 0.0116662 +DEAL:GMRES:: estimated residual at iteration 19: 0.00989870 +DEAL:GMRES:: estimated residual at iteration 20: 0.00800703 +DEAL:GMRES:: estimated residual at iteration 21: 0.00604107 +DEAL:GMRES:: estimated residual at iteration 22: 0.00431736 +DEAL:GMRES:: estimated residual at iteration 23: 0.00304828 +DEAL:GMRES:: estimated residual at iteration 24: 0.00203679 +DEAL:GMRES:: estimated residual at iteration 24: 0.00203679 +DEAL:GMRES:: estimated residual at iteration 25: 0.00141669 +DEAL:GMRES:: estimated residual at iteration 26: 0.00101755 +DEAL:GMRES:: estimated residual at iteration 27: 0.000724887 +DEAL:GMRES:: estimated residual at iteration 28: 0.000535249 +DEAL:GMRES:: estimated residual at iteration 29: 0.000423665 +DEAL:GMRES:: estimated residual at iteration 30: 0.000358414 +DEAL:GMRES:: estimated residual at iteration 30: 0.000358414 +DEAL:GMRES:: estimated residual at iteration 31: 0.000307521 +DEAL:GMRES:: estimated residual at iteration 32: 0.000247031 +DEAL:GMRES:: estimated residual at iteration 33: 0.000183614 +DEAL:GMRES:: estimated residual at iteration 34: 0.000129847 +DEAL:GMRES::Convergence step 35 value 9.20930e-05 +DEAL:GMRES:: estimated residual at iteration 35: 9.20930e-05