From 8ebd39124ec3fec59c50a94b9ec5f88346d2b926 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 6 Jun 2017 23:50:38 +0200 Subject: [PATCH] Use shortcut in multigrid smoothing in case we want to zero out vector --- include/deal.II/multigrid/mg_base.h | 15 ++ include/deal.II/multigrid/mg_smoother.h | 193 ++++++++++++++++-- .../deal.II/multigrid/multigrid.templates.h | 20 +- source/multigrid/mg_base.cc | 11 + 4 files changed, 214 insertions(+), 25 deletions(-) diff --git a/include/deal.II/multigrid/mg_base.h b/include/deal.II/multigrid/mg_base.h index 8e09304e0f..d825ed8b29 100644 --- a/include/deal.II/multigrid/mg_base.h +++ b/include/deal.II/multigrid/mg_base.h @@ -223,6 +223,7 @@ public: * Virtual destructor. */ virtual ~MGSmootherBase(); + /** * Release matrices. */ @@ -234,6 +235,20 @@ public: virtual void smooth (const unsigned int level, VectorType &u, const VectorType &rhs) const = 0; + + /** + * As opposed to the smoothing function, this function applies the action of + * the smoothing, overwriting the previous content in the vector u. This + * function must be equivalent to the following code + * @code + * u = 0; + * smooth(level, u, rhs); + * @endcode + * but can usually be implemented more efficiently than the former. + */ + virtual void apply (const unsigned int level, + VectorType &u, + const VectorType &rhs) const; }; /*@}*/ diff --git a/include/deal.II/multigrid/mg_smoother.h b/include/deal.II/multigrid/mg_smoother.h index 06bbd5632e..11f4236927 100644 --- a/include/deal.II/multigrid/mg_smoother.h +++ b/include/deal.II/multigrid/mg_smoother.h @@ -141,6 +141,7 @@ public: virtual void smooth (const unsigned int level, VectorType &u, const VectorType &rhs) const; + virtual void clear (); }; @@ -224,6 +225,14 @@ namespace mg VectorType &u, const VectorType &rhs) const; + /** + * The apply variant of smoothing, setting the vector u to zero before + * calling the smooth function. + */ + virtual void apply (const unsigned int level, + VectorType &u, + const VectorType &rhs) const; + /** * Memory used by this object. */ @@ -341,6 +350,14 @@ public: VectorType &u, const VectorType &rhs) const; + /** + * The apply variant of smoothing, setting the vector u to zero before + * calling the smooth function. + */ + virtual void apply (const unsigned int level, + VectorType &u, + const VectorType &rhs) const; + /** * Object containing relaxation methods. */ @@ -474,6 +491,14 @@ public: VectorType &u, const VectorType &rhs) const; + /** + * The apply variant of smoothing, setting the vector u to zero before + * calling the smooth function. + */ + virtual void apply (const unsigned int level, + VectorType &u, + const VectorType &rhs) const; + /** * Object containing relaxation methods. */ @@ -654,6 +679,42 @@ namespace mg } + template + inline void + SmootherRelaxation::apply (const unsigned int level, + VectorType &u, + const VectorType &rhs) const + { + unsigned int maxlevel = this->max_level(); + unsigned int steps2 = this->steps; + + if (this->variable) + steps2 *= (1<<(maxlevel-level)); + + bool T = this->transpose; + if (this->symmetric && (steps2 % 2 == 0)) + T = false; + if (this->debug > 0) + deallog << 'S' << level << ' '; + + if (T) + (*this)[level].Tvmult(u, rhs); + else + (*this)[level].vmult(u, rhs); + if (this->symmetric) + T = !T; + for (unsigned int i=1; isymmetric) + T = !T; + } + } + + template inline std::size_t @@ -821,6 +882,42 @@ MGSmootherRelaxation::smooth (const unsi } +template +inline void +MGSmootherRelaxation::apply (const unsigned int level, + VectorType &u, + const VectorType &rhs) const +{ + unsigned int maxlevel = smoothers.max_level(); + unsigned int steps2 = this->steps; + + if (this->variable) + steps2 *= (1<<(maxlevel-level)); + + bool T = this->transpose; + if (this->symmetric && (steps2 % 2 == 0)) + T = false; + if (this->debug > 0) + deallog << 'S' << level << ' '; + + if (T) + smoothers[level].Tvmult(u, rhs); + else + smoothers[level].vmult(u, rhs); + if (this->symmetric) + T = !T; + for (unsigned int i=1; isymmetric) + T = !T; + } +} + + template inline @@ -997,13 +1094,86 @@ MGSmootherPrecondition::smooth { if (this->debug > 0) deallog << 'T'; - if (i == 0 && u.all_zero()) - *r = rhs; - else - { - matrices[level].Tvmult(*r,u); - r->sadd(-1.,1.,rhs); - } + matrices[level].Tvmult(*r,u); + r->sadd(-1.,1.,rhs); + if (this->debug > 2) + deallog << ' ' << r->l2_norm() << ' '; + smoothers[level].Tvmult(*d, *r); + if (this->debug > 1) + deallog << ' ' << d->l2_norm() << ' '; + } + else + { + if (this->debug > 0) + deallog << 'N'; + matrices[level].vmult(*r,u); + r->sadd(-1.,rhs); + if (this->debug > 2) + deallog << ' ' << r->l2_norm() << ' '; + smoothers[level].vmult(*d, *r); + if (this->debug > 1) + deallog << ' ' << d->l2_norm() << ' '; + } + u += *d; + if (this->symmetric) + T = !T; + } + if (this->debug > 0) + deallog << std::endl; +} + + + +template +inline void +MGSmootherPrecondition::apply +(const unsigned int level, + VectorType &u, + const VectorType &rhs) const +{ + unsigned int maxlevel = matrices.max_level(); + unsigned int steps2 = this->steps; + + if (this->variable) + steps2 *= (1<<(maxlevel-level)); + + bool T = this->transpose; + if (this->symmetric && (steps2 % 2 == 0)) + T = false; + if (this->debug > 0) + deallog << 'S' << level << ' '; + + // first step where we overwrite the result + if (this->debug > 2) + deallog << ' ' << rhs.l2_norm() << ' '; + if (this->debug > 0) + deallog << (T ? 'T' : 'N'); + if (T) + smoothers[level].Tvmult(u, rhs); + else + smoothers[level].vmult(u, rhs); + if (this->debug > 1) + deallog << ' ' << u.l2_norm() << ' '; + if (this->symmetric) + T = !T; + + typename VectorMemory::Pointer r(this->vector_memory); + typename VectorMemory::Pointer d(this->vector_memory); + + if (steps2 > 1) + { + r->reinit(u,true); + d->reinit(u,true); + } + + for (unsigned int i=1; idebug > 0) + deallog << 'T'; + matrices[level].Tvmult(*r,u); + r->sadd(-1.,1.,rhs); if (this->debug > 2) deallog << ' ' << r->l2_norm() << ' '; smoothers[level].Tvmult(*d, *r); @@ -1014,13 +1184,8 @@ MGSmootherPrecondition::smooth { if (this->debug > 0) deallog << 'N'; - if (i == 0 && u.all_zero()) - *r = rhs; - else - { - matrices[level].vmult(*r,u); - r->sadd(-1.,rhs); - } + matrices[level].vmult(*r,u); + r->sadd(-1.,rhs); if (this->debug > 2) deallog << ' ' << r->l2_norm() << ' '; smoothers[level].vmult(*d, *r); diff --git a/include/deal.II/multigrid/multigrid.templates.h b/include/deal.II/multigrid/multigrid.templates.h index 833ab779cd..d63f752e49 100644 --- a/include/deal.II/multigrid/multigrid.templates.h +++ b/include/deal.II/multigrid/multigrid.templates.h @@ -125,7 +125,7 @@ Multigrid::level_v_step (const unsigned int level) // smoothing of the residual if (debug>1) deallog << "Smoothing on level " << level << std::endl; - pre_smooth->smooth(level, solution[level], defect[level]); + pre_smooth->apply(level, solution[level], defect[level]); if (debug>2) deallog << "Solution norm " << solution[level].l2_norm() @@ -241,7 +241,7 @@ Multigrid::level_step(const unsigned int level, // smoothing of the residual if (debug>1) deallog << cychar << "-cycle smoothing level " << level << std::endl; - pre_smooth->smooth(level, solution[level], defect2[level]); + pre_smooth->apply(level, solution[level], defect2[level]); if (debug>2) deallog << cychar << "-cycle solution norm " @@ -268,11 +268,7 @@ Multigrid::level_step(const unsigned int level, transfer->restrict_and_add (level, defect2[level-1], t[level]); - // Every cycle starts with a recursion of its type. As opposed to the - // V-cycle above where we do not need to set the solution vector to zero (as - // that has been done in the init functions), we need to zero out the values - // from previous cycles. - solution[level-1] = 0; + // Every cycle starts with a recursion of its type. level_step(level-1, cycle); // For W and F-cycle, repeat the process on the next coarser level except @@ -337,8 +333,10 @@ Multigrid::cycle() for (unsigned int level=minlevel; level<=maxlevel; ++level) { - solution[level].reinit(defect[level]); - t[level].reinit(defect[level]); + // the vectors for level>minlevel will be overwritten by the apply() + // method of the smoother -> do not force them to be zeroed out here + solution[level].reinit(defect[level], level>minlevel); + t[level].reinit(defect[level], level>minlevel); if (cycle_type != v_cycle) defect2[level].reinit(defect[level]); } @@ -362,8 +360,8 @@ Multigrid::vcycle() for (unsigned int level=minlevel; level<=maxlevel; ++level) { - solution[level].reinit(defect[level]); - t[level].reinit(defect[level]); + solution[level].reinit(defect[level], level>minlevel); + t[level].reinit(defect[level], level>minlevel); } level_v_step (maxlevel); } diff --git a/source/multigrid/mg_base.cc b/source/multigrid/mg_base.cc index 37fa24b6b6..efb63321dd 100644 --- a/source/multigrid/mg_base.cc +++ b/source/multigrid/mg_base.cc @@ -42,6 +42,17 @@ MGSmootherBase::~MGSmootherBase() {} +template +void +MGSmootherBase::apply (const unsigned int level, + VectorType &u, + const VectorType &rhs) const +{ + u = 0; + smooth(level, u, rhs); +} + + template MGCoarseGridBase::~MGCoarseGridBase() {} -- 2.39.5