From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Tue, 6 Jun 2017 21:50:38 +0000 (+0200)
Subject: Use shortcut in multigrid smoothing in case we want to zero out vector
X-Git-Tag: v9.0.0-rc1~1513^2~3
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8ebd39124ec3fec59c50a94b9ec5f88346d2b926;p=dealii.git

Use shortcut in multigrid smoothing in case we want to zero out vector
---

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 <class RelaxationType, typename VectorType>
+  inline void
+  SmootherRelaxation<RelaxationType, VectorType>::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; i<steps2; ++i)
+      {
+        if (T)
+          (*this)[level].Tstep(u, rhs);
+        else
+          (*this)[level].step(u, rhs);
+        if (this->symmetric)
+          T = !T;
+      }
+  }
+
+
   template <class RelaxationType, typename VectorType>
   inline
   std::size_t
@@ -821,6 +882,42 @@ MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::smooth (const unsi
 }
 
 
+template <typename MatrixType, class RelaxationType, typename VectorType>
+inline void
+MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::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; i<steps2; ++i)
+    {
+      if (T)
+        smoothers[level].Tstep(u, rhs);
+      else
+        smoothers[level].step(u, rhs);
+      if (this->symmetric)
+        T = !T;
+    }
+}
+
+
 
 template <typename MatrixType, class RelaxationType, typename VectorType>
 inline
@@ -997,13 +1094,86 @@ MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::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 <typename MatrixType, typename PreconditionerType, typename VectorType>
+inline void
+MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::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<VectorType>::Pointer r(this->vector_memory);
+  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
+
+  if (steps2 > 1)
+    {
+      r->reinit(u,true);
+      d->reinit(u,true);
+    }
+
+  for (unsigned int i=1; i<steps2; ++i)
+    {
+      if (T)
+        {
+          if (this->debug > 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<MatrixType, PreconditionerType, VectorType>::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<VectorType>::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<VectorType>::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<VectorType>::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<VectorType>::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<VectorType>::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<VectorType>::~MGSmootherBase()
 {}
 
 
+template <typename VectorType>
+void
+MGSmootherBase<VectorType>::apply (const unsigned int level,
+                                   VectorType         &u,
+                                   const VectorType   &rhs) const
+{
+  u = 0;
+  smooth(level, u, rhs);
+}
+
+
 template <typename VectorType>
 MGCoarseGridBase<VectorType>::~MGCoarseGridBase()
 {}