From: Martin Kronbichler Date: Fri, 21 Oct 2016 14:18:23 +0000 (+0200) Subject: Add step() and Tstep() methods to Chebyshev preconditioner. X-Git-Tag: v8.5.0-rc1~547^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=34e48106cfd096ef05c35df2edc2eea19df969e3;p=dealii.git Add step() and Tstep() methods to Chebyshev preconditioner. --- diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 6998e0f332..dface084ed 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -877,6 +877,9 @@ public: * smoothing in a multigrid algorithm), but not the way the solver classes * expect a preconditioner to work (where one ignores the content in * dst for the preconditioner application). + * + * @deprecated For non-zero starting, use the step() and Tstep() + * interfaces, whereas vmult() provides the preconditioner interface. */ bool nonzero_starting; @@ -936,6 +939,18 @@ public: void Tvmult (VectorType &dst, const VectorType &src) const; + /** + * Perform one step of the preconditioned Richardson iteration. + */ + void step (VectorType &dst, + const VectorType &src) const; + + /** + * Perform one transposed step of the preconditioned Richardson iteration. + */ + void Tstep (VectorType &dst, + const VectorType &src) const; + /** * Resets the preconditioner. */ @@ -990,6 +1005,22 @@ private: * Stores whether the preconditioner has been set up. */ bool is_initialized; + + /** + * Runs the inner loop of the Chebyshev preconditioner that is the same for + * vmult() and step() methods. + */ + void do_chebyshev_loop(VectorType &dst, + const VectorType &src) const; + + /** + * Runs the inner loop of the Chebyshev preconditioner that is the same for + * vmult() and step() methods. Uses a separate function to not force users + * to provide both vmult() and Tvmult() in case only one variant is + * requested in subsequent calls. + */ + void do_transpose_chebyshev_loop(VectorType &dst, + const VectorType &src) const; }; @@ -1913,26 +1944,34 @@ PreconditionChebyshev::initialize template inline void -PreconditionChebyshev::vmult (VectorType &dst, - const VectorType &src) const +PreconditionChebyshev::do_chebyshev_loop(VectorType &dst, + const VectorType &src) const { - Assert (is_initialized, ExcMessage("Preconditioner not initialized")); double rhok = delta / theta, sigma = theta / delta; - if (data.nonzero_starting && !dst.all_zero()) + for (unsigned int k=0; kvmult (update2, dst); + const double rhokp = 1./(2.*sigma-rhok); + const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta; + rhok = rhokp; internal::PreconditionChebyshev::vector_updates - (src, data.matrix_diagonal_inverse, false, 0., 1./theta, update1, + (src, data.matrix_diagonal_inverse, false, factor1, factor2, update1, update2, dst); } - else - internal::PreconditionChebyshev::vector_updates - (src, data.matrix_diagonal_inverse, true, 0., 1./theta, update1, - update2, dst); +} + + +template +inline +void +PreconditionChebyshev::do_transpose_chebyshev_loop(VectorType &dst, + const VectorType &src) const +{ + double rhok = delta / theta, sigma = theta / delta; for (unsigned int k=0; kvmult (update2, dst); + matrix_ptr->Tvmult (update2, dst); const double rhokp = 1./(2.*sigma-rhok); const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta; rhok = rhokp; @@ -1944,6 +1983,22 @@ PreconditionChebyshev::vmult (VectorType &dst, +template +inline +void +PreconditionChebyshev::vmult (VectorType &dst, + const VectorType &src) const +{ + Assert (is_initialized, ExcMessage("Preconditioner not initialized")); + internal::PreconditionChebyshev::vector_updates + (src, data.matrix_diagonal_inverse, true, 0., 1./theta, update1, + update2, dst); + + do_chebyshev_loop(dst, src); +} + + + template inline void @@ -1951,10 +2006,25 @@ PreconditionChebyshev::Tvmult (VectorType &dst, const VectorType &src) const { Assert (is_initialized, ExcMessage("Preconditioner not initialized")); - double rhok = delta / theta, sigma = theta / delta; - if (data.nonzero_starting && !dst.all_zero()) + internal::PreconditionChebyshev::vector_updates + (src, data.matrix_diagonal_inverse, true, 0., 1./theta, update1, + update2, dst); + + do_transpose_chebyshev_loop(dst, src); +} + + + +template +inline +void +PreconditionChebyshev::step (VectorType &dst, + const VectorType &src) const +{ + Assert (is_initialized, ExcMessage("Preconditioner not initialized")); + if (!dst.all_zero()) { - matrix_ptr->Tvmult (update2, dst); + matrix_ptr->vmult (update2, dst); internal::PreconditionChebyshev::vector_updates (src, data.matrix_diagonal_inverse, false, 0., 1./theta, update1, update2, dst); @@ -1964,16 +2034,31 @@ PreconditionChebyshev::Tvmult (VectorType &dst, (src, data.matrix_diagonal_inverse, true, 0., 1./theta, update1, update2, dst); - for (unsigned int k=0; k +inline +void +PreconditionChebyshev::Tstep (VectorType &dst, + const VectorType &src) const +{ + Assert (is_initialized, ExcMessage("Preconditioner not initialized")); + if (!dst.all_zero()) { matrix_ptr->Tvmult (update2, dst); - const double rhokp = 1./(2.*sigma-rhok); - const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta; - rhok = rhokp; internal::PreconditionChebyshev::vector_updates - (src, data.matrix_diagonal_inverse, false, factor1, factor2, update1, + (src, data.matrix_diagonal_inverse, false, 0., 1./theta, update1, update2, dst); } + else + internal::PreconditionChebyshev::vector_updates + (src, data.matrix_diagonal_inverse, true, 0., 1./theta, update1, + update2, dst); + + do_transpose_chebyshev_loop(dst, src); }