trying to find the bottlenec in copy_to_mg

git-svn-id: https://svn.dealii.org/trunk@9586 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/deal.II/include/multigrid/mg_transfer.templates.h b/deal.II/deal.II/include/multigrid/mg_transfer.templates.h
index 54e53b23c35515dc787d55c2312a43019c2caa52..15301e7f752750b2dfa91943450098e1d62ef41b 100644 (file)
--- a/deal.II/deal.II/include/multigrid/mg_transfer.templates.h
+++ b/deal.II/deal.II/include/multigrid/mg_transfer.templates.h
@@ -362,6 +362,18 @@ MGTransferSelect<number>::do_copy_to_mg (
   std::vector<unsigned int> level_dof_indices  (dofs_per_cell);
   std::vector<unsigned int> level_face_indices (dofs_per_face);
 
+                                  // Build a vector of the selected
+                                  // indices, since traversing all
+                                  // indices on each cell is too
+                                  // slow.
+  std::vector<unsigned int> selected_indices;
+  selected_indices.reserve(dofs_per_cell);
+  for (unsigned int i=0; i<dofs_per_cell; ++i)
+    if (mg_selected[mg_target_component[fe.system_to_component_index(i).first]])
+      selected_indices.push_back(i);
+  unsigned int selected_component
+    = mg_target_component[fe.system_to_component_index(selected_indices[0]).first];
+
                                   // traverse the grid top-down
                                   // (i.e. starting with the most
                                   // refined grid). this way, we can
@@ -398,36 +410,18 @@ MGTransferSelect<number>::do_copy_to_mg (
                                           // transfer the global
                                           // defect in the vector
                                           // into the level-wise one
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
+         const unsigned int level_start
+           = mg_component_start[level][selected_component];
+         const typename std::vector<unsigned int>::const_iterator
+           end = selected_indices.end();
+         
+         for (typename std::vector<unsigned int>::const_iterator
+                i=selected_indices.begin();
+              i != end; ++i)
            {
-             const unsigned int component
-               = mg_target_component[fe.system_to_component_index(i).first];
-             if (mg_selected[component])
-               {
-                 const unsigned int level_start
-                   = mg_component_start[level][component];
-                 
-                 dst[level](level_dof_indices[i] - level_start)
-                   = src(global_dof_indices[i] - offset);
-               }
+             dst[level](level_dof_indices[*i] - level_start)
+               = src(global_dof_indices[*i] - offset);
            }
-         
-//       for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
-//         {
-//           const typename MGDoFHandler<dim>::face_iterator
-//             face = level_cell->face(face_n);
-//           if (face->has_children())
-//             {
-//               face->get_mg_dof_indices(level_face_indices);
-
-
-                                                  // Delete values on refinement edge,
-                                                  // since restriction will add them again.
-//               for (unsigned int i=0; i<dofs_per_face; ++i)
-//                 dst[level](level_face_indices[i])
-//                   = 0.;
-//             };
-//         };
        }
                                       // for that part of the level
                                       // which is further refined:
@@ -436,7 +430,7 @@ MGTransferSelect<number>::do_copy_to_mg (
                                       // one level higher
       if (static_cast<unsigned int>(level) < maxlevel)
        {
-         restrict_and_add (level+1, dst[level], dst[level+1]);
+         ;//     restrict_and_add (level+1, dst[level], dst[level+1]);
        }
     };
 }

diff --git a/deal.II/deal.II/include/multigrid/multigrid.h b/deal.II/deal.II/include/multigrid/multigrid.h
index 2510b504f4ad6810552d3b0a9a128e62d43bfba3..94ee218b8bd21dc1b115314ae0c1d8e3f148fdc9 100644 (file)
--- a/deal.II/deal.II/include/multigrid/multigrid.h
+++ b/deal.II/deal.II/include/multigrid/multigrid.h
@@ -402,12 +402,18 @@ PreconditionMG<dim, VECTOR, TRANSFER>::vmult (
   VECTOR2& dst,
   const VECTOR2& src) const
 {
+  deallog << "copy to" << endl;
+  
   transfer->copy_to_mg(*mg_dof_handler,
                       multigrid->defect,
                       src);
   
+  deallog << "mg" << endl;
+  
   multigrid->vcycle();
 
+  deallog << "copy from " << endl;
+  
   transfer->copy_from_mg(*mg_dof_handler,
                         dst,
                         multigrid->solution);
@@ -421,12 +427,18 @@ PreconditionMG<dim, VECTOR, TRANSFER>::vmult_add (
   VECTOR2& dst,
   const VECTOR2& src) const
 {
+  deallog << "copy to" << endl;
+  
   transfer->copy_to_mg(*mg_dof_handler,
                       multigrid->defect,
                       src);
+
+  deallog << "mg" << endl;
   
   multigrid->vcycle();
 
+  deallog << "copy from " << endl;
+  
   transfer->copy_from_mg_add(*mg_dof_handler,
                             dst,
                             multigrid->solution);

diff --git a/deal.II/deal.II/include/multigrid/multigrid.templates.h b/deal.II/deal.II/include/multigrid/multigrid.templates.h
index 53fd6f42d5849582f80e2a41a05c0a724708be4c..183f6df8448203540389a01c3c2e1326fb5e13dc 100644 (file)
--- a/deal.II/deal.II/include/multigrid/multigrid.templates.h
+++ b/deal.II/deal.II/include/multigrid/multigrid.templates.h
@@ -14,6 +14,11 @@
 #define __deal2__multigrid_templates_h
 #include <multigrid/multigrid.h>
 
+#include <base/logstream.h>
+
+#include "iostream"
+
+using namespace std;
 
 template <class VECTOR>
 void
@@ -39,13 +44,17 @@ Multigrid<VECTOR>::level_mgstep(const unsigned int level)
       return;
     }
 
+//  deallog << "Pre-smooth " << level << endl;
+  
                                   // smoothing of the residual by
                                   // modifying s
   pre_smooth->smooth(level, solution[level], defect[level]);
 
+//  deallog << "vmult " << level << endl;
                                   // t = A*solution[level]
   matrix->vmult(level, t[level], solution[level]);
-  
+
+//  deallog << "restrict " << level << endl;
                                   // make t rhs of lower level The
                                   // non-refined parts of the
                                   // coarse-level defect already
@@ -62,6 +71,7 @@ Multigrid<VECTOR>::level_mgstep(const unsigned int level)
       defect[l-1] -= t[l-1];
     }
 
+//    deallog << "recursion " << level << endl;
                                   // do recursion
   level_mgstep(level-1);
 
@@ -72,7 +82,7 @@ Multigrid<VECTOR>::level_mgstep(const unsigned int level)
   t[level] = 0.;
 
                                   // do coarse grid correction
-
+//  deallog << "prolongate " << level << endl;
   transfer->prolongate(level, t[level], solution[level-1]);
   
   solution[level] += t[level];
@@ -82,9 +92,10 @@ Multigrid<VECTOR>::level_mgstep(const unsigned int level)
       edge_up->Tvmult(level, t[level], solution[level-1]);
       defect[level] -= t[level];
     }
-  
+//  deallog << "Post-smooth " << level << endl;
                                   // post-smoothing
   post_smooth->smooth(level, solution[level], defect[level]);
+//  deallog << "ready " << level << endl;
 }
 
 

diff --git a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc
index fab658d581ca144aaef53c34a477bd2eb2a88004..d8d4bbf49df0523e7c7e68d0f7eaf911f9657108 100644 (file)
--- a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc
+++ b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -23,7 +23,7 @@
 #include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
 
-
+//TODO:[GK] Inline simple functions in 1d and 3d
 
 /* ------------------------ MGDoFLineAccessor --------------------------- */
 
@@ -210,6 +210,7 @@ MGDoFObjectAccessor<2, dim>::MGDoFObjectAccessor (const Triangulation<dim> *tria
 
 
 template <int dim>
+inline
 unsigned int MGDoFObjectAccessor<2, dim>::mg_dof_index (const unsigned int i) const
 {
   Assert (this->dof_handler != 0,
@@ -229,6 +230,7 @@ unsigned int MGDoFObjectAccessor<2, dim>::mg_dof_index (const unsigned int i) co
 
 
 template <int dim>
+inline
 void MGDoFObjectAccessor<2, dim>::set_mg_dof_index (const unsigned int i,
                                                    const unsigned int index) const
 {
@@ -249,6 +251,7 @@ void MGDoFObjectAccessor<2, dim>::set_mg_dof_index (const unsigned int i,
 
 
 template <int dim>
+inline
 unsigned int MGDoFObjectAccessor<2, dim>::mg_vertex_dof_index (const unsigned int vertex,
                                                               const unsigned int i) const
 {
@@ -268,6 +271,7 @@ unsigned int MGDoFObjectAccessor<2, dim>::mg_vertex_dof_index (const unsigned in
 
 
 template <int dim>
+inline
 void MGDoFObjectAccessor<2, dim>::set_mg_vertex_dof_index (const unsigned int vertex,
                                                           const unsigned int i,
                                                           const unsigned int index) const
@@ -353,6 +357,7 @@ MGDoFObjectAccessor<2,dim>::get_mg_dof_values (const Vector<number> &values,
 
 
 template <int dim>
+inline
 TriaIterator<dim,MGDoFObjectAccessor<1, dim> >
 MGDoFObjectAccessor<2, dim>::line (const unsigned int i) const
 {
@@ -369,6 +374,7 @@ MGDoFObjectAccessor<2, dim>::line (const unsigned int i) const
 
 
 template <int dim>
+inline
 TriaIterator<dim,MGDoFObjectAccessor<2, dim> >
 MGDoFObjectAccessor<2, dim>::child (const unsigned int i) const
 {

diff --git a/deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc b/deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc
index 18528d514fe6c09bb088ddd6ab4ad94e91b3be11..e6d9ce714d1f84f73d4fde3edf41a546f20d9a03 100644 (file)
--- a/deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc
+++ b/deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc
@@ -17,6 +17,7 @@
 #include <lac/sparse_matrix.h>
 #include <lac/block_sparse_matrix.h>
 #include <multigrid/mg_transfer.h>
+#include <multigrid/mg_transfer.templates.h>
 #include <multigrid/multigrid.templates.h>
 
 

diff --git a/deal.II/lac/source/sparsity_pattern.cc b/deal.II/lac/source/sparsity_pattern.cc
index 37d4503742b60fd1f7035d8db34cd11a8b059bbe..96657d9ddd892e00610472d3d33418265937e79f 100644 (file)
--- a/deal.II/lac/source/sparsity_pattern.cc
+++ b/deal.II/lac/source/sparsity_pattern.cc
@@ -330,7 +330,12 @@ SparsityPattern::reinit (const unsigned int               m,
                                   // exactly one element, to have a
                                   // valid pointer to some memory
   if (vec_len == 0)
-    vec_len = 1;
+    {
+      vec_len = 1;
+      if (colnums) delete[] colnums;
+      max_vec_len = vec_len;
+      colnums = new unsigned int[max_vec_len];
+    }
   
   max_row_length = (row_lengths.size() == 0 ?
                    0 :
@@ -339,7 +344,7 @@ SparsityPattern::reinit (const unsigned int               m,
 
   if (diagonal_optimized && (max_row_length==0) && (m!=0))
     max_row_length = 1;
-
+  
                                   // allocate memory for the rowstart
                                   // values, if necessary
   if (rows > max_dim)