From: Martin Kronbichler Date: Wed, 11 Jan 2017 14:58:28 +0000 (+0100) Subject: Use common vector interface for Chebyshev smoother X-Git-Tag: v8.5.0-rc1~269^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3e8227c930a50d63c377eac164784b2a561ac52d;p=dealii.git Use common vector interface for Chebyshev smoother --- diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 6af9b36195..c201cf15cd 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -1903,7 +1903,6 @@ namespace internal (void)matrix; (void)preconditioner; AssertThrow(preconditioner.get() != NULL, ExcNotInitialized()); - AssertDimension(preconditioner->m(), matrix.m()); } template @@ -1925,7 +1924,10 @@ namespace internal // Check if we can initialize from vector that then gets set to zero // as the matrix will own the memory preconditioner->reinit(diagonal_inverse); - diagonal_inverse.reinit(0); + { + VectorType empty_vector; + diagonal_inverse.reinit(empty_vector); + } // This part only works in serial if (preconditioner->m() != matrix.m()) @@ -2050,11 +2052,14 @@ PreconditionChebyshev::clear () is_initialized = false; theta = delta = 1.0; matrix_ptr = 0; - data.matrix_diagonal_inverse.reinit(0); + { + VectorType empty_vector; + data.matrix_diagonal_inverse.reinit(empty_vector); + update1.reinit(empty_vector); + update2.reinit(empty_vector); + update3.reinit(empty_vector); + } data.preconditioner.reset(); - update1.reinit(0); - update2.reinit(0); - update3.reinit(0); } DEAL_II_ENABLE_EXTRA_DIAGNOSTICS diff --git a/tests/lac/precondition_chebyshev_05.cc b/tests/lac/precondition_chebyshev_05.cc new file mode 100644 index 0000000000..97ba24a229 --- /dev/null +++ b/tests/lac/precondition_chebyshev_05.cc @@ -0,0 +1,100 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test PreconditionChebyshev on parallel Trilinos matrix and vector + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include "../testmatrix.h" +#include + +#include +#include +#include +#include + + + +int main(int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog << std::fixed; + deallog << std::setprecision(4); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + + + for (unsigned int size=4; size <= 16; size *= 2) + { + unsigned int dim = (size-1)*(size-1); + + deallog << "Size " << size << " Unknowns " << dim << std::endl; + + // Make matrix + FDMatrix testproblem(size, size); + SparsityPattern structure(dim, dim, 5); + testproblem.five_point_structure(structure); + structure.compress(); + SparseMatrix A(structure); + testproblem.five_point(A); + TrilinosWrappers::SparseMatrix AA; + AA.reinit(A); + + PreconditionChebyshev cheby; + PreconditionChebyshev::AdditionalData cheby_data; + cheby_data.preconditioner.reset(new TrilinosWrappers::PreconditionJacobi()); + cheby_data.preconditioner->initialize(AA); + cheby_data.degree = 10; + cheby_data.smoothing_range = 40; + cheby.initialize(AA, cheby_data); + + IndexSet set(dim); + set.add_range(0, dim); + TrilinosWrappers::MPI::Vector v, tmp1, tmp2; + v.reinit(set, MPI_COMM_WORLD); + tmp1.reinit(set, MPI_COMM_WORLD); + tmp2.reinit(set, MPI_COMM_WORLD); + for (unsigned int i=0; i<3; ++i) + { + for (unsigned int j=0; jvmult (tmp2, tmp1); + tmp2 -= v; + const double ilu_residual = tmp2.l2_norm(); + + AA.vmult (tmp1, v); + cheby.vmult (tmp2, tmp1); + tmp2 -= v; + const double cheby_residual = tmp2.l2_norm(); + + deallog << "Residual step i=" << i << ": " + << " jacobi=" << ilu_residual + << ", cheby=" << cheby_residual + << std::endl; + } + } + + return 0; +} diff --git a/tests/lac/precondition_chebyshev_05.with_trilinos=true.with_lapack=true.output b/tests/lac/precondition_chebyshev_05.with_trilinos=true.with_lapack=true.output new file mode 100644 index 0000000000..7102557826 --- /dev/null +++ b/tests/lac/precondition_chebyshev_05.with_trilinos=true.with_lapack=true.output @@ -0,0 +1,13 @@ + +DEAL::Size 4 Unknowns 9 +DEAL::Residual step i=0: jacobi=1.2946, cheby=0.0432 +DEAL::Residual step i=1: jacobi=1.3005, cheby=0.0339 +DEAL::Residual step i=2: jacobi=0.7640, cheby=0.0260 +DEAL::Size 8 Unknowns 49 +DEAL::Residual step i=0: jacobi=3.3132, cheby=0.1465 +DEAL::Residual step i=1: jacobi=3.5095, cheby=0.1507 +DEAL::Residual step i=2: jacobi=2.9701, cheby=0.1604 +DEAL::Size 16 Unknowns 225 +DEAL::Residual step i=0: jacobi=7.5744, cheby=3.4216 +DEAL::Residual step i=1: jacobi=7.6760, cheby=3.4354 +DEAL::Residual step i=2: jacobi=7.0389, cheby=3.1193