]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Introduce a field that stores which constraints exist already. This is desirable...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 4 Oct 2004 21:29:55 +0000 (21:29 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 4 Oct 2004 21:29:55 +0000 (21:29 +0000)
git-svn-id: https://svn.dealii.org/trunk@9689 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/source/dofs/dof_constraints.cc
deal.II/doc/news/c-5.0.html

index 8e46abc76324f01d4c0a1e560c5db132cd5aa2d5..77c55610455f9b30e7222cd17491d56536cc24dc 100644 (file)
@@ -898,6 +898,53 @@ class ConstraintMatrix : public Subscriptor
                                      * much slower.
                                      */
     std::vector<ConstraintLine> lines;
+
+                                    /**
+                                     * A list of flags that indicate
+                                     * whether there is a constraint
+                                     * line for a given degree of
+                                     * freedom index. Note that this
+                                     * class has no notion of how
+                                     * many degrees of freedom there
+                                     * really are, so if we check
+                                     * whether there is a constraint
+                                     * line for a given degree of
+                                     * freedom, then this vector may
+                                     * actually be shorter than the
+                                     * index of the DoF we check for.
+                                     *
+                                     * This field exists since when adding a
+                                     * new constraint line we have to figure
+                                     * out whether it already
+                                     * exists. Previously, we would simply
+                                     * walk the unsorted list of constraint
+                                     * lines until we either hit the end or
+                                     * found it. This algorithm is O(N) if N
+                                     * is the number of constraints, which
+                                     * makes it O(N^2) when inserting all
+                                     * constraints. For large problems with
+                                     * many constraints, this could easily
+                                     * take 5-10 per cent of the total run
+                                     * time. With this field, we can at least
+                                     * save this time when checking whether a
+                                     * new constraint line already exists.
+                                     *
+                                     * To make things worse, traversing the
+                                     * list of existing constraints requires
+                                     * reads from many different places in
+                                     * memory. Thus, in large 3d
+                                     * applications, the add_line() function
+                                     * showed up very prominently in the
+                                     * overall compute time, mainly because
+                                     * it generated a lot of cache
+                                     * misses. This should also be fixed by
+                                     * using the O(1) algorithm to access the
+                                     * fields of this array.
+                                     * 
+                                     * The field is useful in a number of
+                                     * other contexts as well, though.
+                                     */
+    std::vector<bool> constraint_line_exists;
        
                                     /**
                                      * Store whether the arrays are sorted.
@@ -917,4 +964,101 @@ class ConstraintMatrix : public Subscriptor
 };
 
 
+
+/* ---------------- template and inline functions ----------------- */
+
+inline
+void
+ConstraintMatrix::add_line (const unsigned int line)
+{
+  Assert (sorted==false, ExcMatrixIsClosed());
+
+  
+  
+                                  // check whether line already
+                                  // exists; it may, but then we need
+                                  // can just quit
+  if ((line < constraint_line_exists.size())
+      &&
+      (constraint_line_exists[line] == true))
+    return;
+
+                                  // if necessary enlarge vector of
+                                  // existing entries
+  if (line >= constraint_line_exists.size())
+    constraint_line_exists.resize (line+1, false);
+  constraint_line_exists[line] = true;
+  
+                                  // push a new line to the end of the
+                                  // list
+  lines.push_back (ConstraintLine());
+  lines.back().line = line;
+}
+
+
+
+inline
+void
+ConstraintMatrix::add_entry (const unsigned int line,
+                             const unsigned int column,
+                             const double       value)
+{
+  Assert (sorted==false, ExcMatrixIsClosed());
+
+  std::vector<ConstraintLine>::iterator line_ptr;
+  const std::vector<ConstraintLine>::const_iterator start=lines.begin();
+
+                                  // the usual case is that the line where
+                                  // a value is entered is the one we
+                                  // added last, so we search backward
+  for (line_ptr=(lines.end()-1); line_ptr!=start; --line_ptr)
+    if (line_ptr->line == line)
+      break;
+
+                                  // if the loop didn't break, then
+                                  // line_ptr must be begin().
+                                  // we have an error if that doesn't
+                                  // point to 'line' then
+  Assert (line_ptr->line==line, ExcLineInexistant(line));
+                                  // second check: the respective
+                                  // flag in the
+                                  // constraint_line_exists field
+                                  // must exist
+  Assert (line < constraint_line_exists.size(), ExcInternalError());
+  Assert (constraint_line_exists[line] == true, ExcInternalError());
+  
+
+                                  // if in debug mode, check whether an
+                                  // entry for this column already
+                                  // exists and if its the same as
+                                  // the one entered at present
+                                  //
+                                  // in any case: exit the function if an
+                                  // entry for this column already exists,
+                                  // since we don't want to enter it twice
+  for (std::vector<std::pair<unsigned int,double> >::const_iterator
+         p=line_ptr->entries.begin();
+       p != line_ptr->entries.end(); ++p)
+    if (p->first == column)
+      {
+       Assert (p->second == value,
+               ExcEntryAlreadyExists(line, column, p->second, value));
+       return;
+      }
+  
+  line_ptr->entries.push_back (std::make_pair(column,value));
+}
+
+
+
+inline
+bool
+ConstraintMatrix::is_constrained (const unsigned int index) const 
+{
+  return ((index < constraint_line_exists.size())
+          &&
+          (constraint_line_exists[index] == true));
+}
+
+
 #endif
index d35ca48ff7b3a1bedd18236db311db317ed037fe..51d0598b0b3792f4d95c0f673b71796fa99f4daf 100644 (file)
@@ -80,75 +80,13 @@ ConstraintMatrix::ConstraintLine::memory_consumption () const
 
 
 
-ConstraintMatrix::ConstraintMatrix () :
-               lines(),
-               sorted(false)
+ConstraintMatrix::ConstraintMatrix ()
+               :
+               lines (),
+               sorted (false)
 {}
 
 
-void ConstraintMatrix::add_line (const unsigned int line)
-{
-  Assert (sorted==false, ExcMatrixIsClosed());
-
-                                  // check whether line already exists;
-                                  // it may, but then we need to quit
-  for (unsigned int i=0; i!=lines.size(); ++i)
-    if (lines[i].line == line)
-      return;
-
-                                  // push a new line to the end of the
-                                  // list
-  lines.push_back (ConstraintLine());
-  lines.back().line = line;
-}
-
-
-
-void
-ConstraintMatrix::add_entry (const unsigned int line,
-                             const unsigned int column,
-                             const double       value)
-{
-  Assert (sorted==false, ExcMatrixIsClosed());
-
-  std::vector<ConstraintLine>::iterator line_ptr;
-  const std::vector<ConstraintLine>::const_iterator start=lines.begin();
-                                  // the usual case is that the line where
-                                  // a value is entered is the one we
-                                  // added last, so we search backward
-  for (line_ptr=(lines.end()-1); line_ptr!=start; --line_ptr)
-    if (line_ptr->line == line)
-      break;
-
-                                  // if the loop didn't break, then
-                                  // line_ptr must be begin().
-                                  // we have an error if that doesn't
-                                  // point to 'line' then
-  Assert (line_ptr->line==line, ExcLineInexistant(line));
-
-                                  // if in debug mode, check whether an
-                                  // entry for this column already
-                                  // exists and if its the same as
-                                  // the one entered at present
-                                  //
-                                  // in any case: exit the function if an
-                                  // entry for this column already exists,
-                                  // since we don't want to enter it twice
-  for (std::vector<std::pair<unsigned int,double> >::const_iterator
-         p=line_ptr->entries.begin();
-       p != line_ptr->entries.end(); ++p)
-    if (p->first == column)
-      {
-       Assert (p->second == value,
-               ExcEntryAlreadyExists(line, column, p->second, value));
-       return;
-      };
-  
-  line_ptr->entries.push_back (std::make_pair(column,value));
-}
-
-
-
 void
 ConstraintMatrix::add_entries (const unsigned int                        line,
                                const std::vector<std::pair<unsigned int,double> > &col_val_pairs)
@@ -341,6 +279,15 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints)
   const bool object_was_sorted = sorted;
   sorted = false;
 
+                                  // before we even start: merge the
+                                  // two flag arrays
+  if (other_constraints.constraint_line_exists.size() >
+      constraint_line_exists.size())
+    constraint_line_exists.resize (other_constraints.constraint_line_exists.size(),
+                                  false);
+  for (unsigned int i=0; i<other_constraints.constraint_line_exists.size(); ++i)
+    if (other_constraints.constraint_line_exists[i] == true)
+      constraint_line_exists[i] = true;
 
                                   // first action is to fold into the
                                   // present object possible
@@ -499,6 +446,9 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints)
 
 void ConstraintMatrix::shift (const unsigned int offset)
 {
+  constraint_line_exists.insert (constraint_line_exists.begin(), offset,
+                                false);
+  
   for (std::vector<ConstraintLine>::iterator i = lines.begin();
        i != lines.end(); i++)
     {
@@ -514,8 +464,16 @@ void ConstraintMatrix::shift (const unsigned int offset)
 
 void ConstraintMatrix::clear ()
 {
-  std::vector<ConstraintLine> tmp;
-  lines.swap (tmp);
+  {
+    std::vector<ConstraintLine> tmp;
+    lines.swap (tmp);
+  }
+  
+  {
+    std::vector<bool> tmp;
+    constraint_line_exists.swap (tmp);
+  }
+  
   sorted = false;
 }
 
@@ -1139,36 +1097,9 @@ unsigned int ConstraintMatrix::n_constraints () const
 
 
 
-bool ConstraintMatrix::is_constrained (const unsigned int index) const 
-{
-  if (lines.size() == 0)
-    return false;
-  
-  if (sorted == true)
-    {
-      ConstraintLine index_comparison;
-      index_comparison.line = index;
-      
-      return std::binary_search (lines.begin (),
-                                lines.end (),
-                                index_comparison);
-    }
-  else
-    {
-      for (std::vector<ConstraintLine>::const_iterator i=lines.begin();
-          i!=lines.end(); ++i)
-       if (i->line == index)
-         return true;
-
-      return false;
-    };
-}
-
-
-
 bool ConstraintMatrix::is_identity_constrained (const unsigned int index) const 
 {
-  if (lines.size() == 0)
+  if (is_constrained(index) == false)
     return false;
   
   if (sorted == true)
@@ -1241,6 +1172,7 @@ unsigned int
 ConstraintMatrix::memory_consumption () const
 {
   return (MemoryConsumption::memory_consumption (lines) +
+         MemoryConsumption::memory_consumption (constraint_line_exists) +
          MemoryConsumption::memory_consumption (sorted));
 }
 
index 9a1094dac1e56aa638a5d3f6ce4b24b5d19ace29..81a95a6bf99af1b082584bf97a798426476322ca 100644 (file)
@@ -343,6 +343,19 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       Fixed: The <code class="class">ConstraintMatrix</code> class had some
+       algorithms that were linear in the number of constraints. Since these
+       functions had to be called for each constraint, this resulted in a
+       quadratic behavior. To make things worse, these algorithms traversed
+       large memory blocks leading to a vast number of cache misses which made
+       them really slow. This is now fixed: the algorithm is O(1) and should
+       only access single elements in memory, and one 3d testcase is now a
+       full 5 per cent faster on about 10 minutes runtime.
+       <br>
+       (WB 2004/10/04)
+       </p>
+
   <li> <p>
        New: The <code class="member">GridOut::write_gnuplot</code>
        function now supports curved boundaries also for <code>dim==3</code>.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.