From 34a881bbeecf5b3fcca7fbba18433dd7448f6641 Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 21 Jun 2004 17:02:47 +0000
Subject: [PATCH] Use a more sophisticated algorithm for the compressed
 sparsity pattern that uses caches and std::vector instead of std::set.
 std::set allocates 20 bytes each time we add a single integer, which can't be
 efficient both from memory and run-time viewpoints.

git-svn-id: https://svn.dealii.org/trunk@9440 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../include/lac/compressed_sparsity_pattern.h | 306 +++++++++++++++++-
 .../lac/source/compressed_sparsity_pattern.cc | 120 ++++---
 2 files changed, 362 insertions(+), 64 deletions(-)

diff --git a/deal.II/lac/include/lac/compressed_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_sparsity_pattern.h
index 4e675de5a5..3ad0248d9e 100644
--- a/deal.II/lac/include/lac/compressed_sparsity_pattern.h
+++ b/deal.II/lac/include/lac/compressed_sparsity_pattern.h
@@ -347,12 +347,93 @@ class CompressedSparsityPattern : public Subscriptor
 				      */
     unsigned int cols;
 
+                                     /**
+                                      * Store some data for each row
+                                      * describing which entries of this row
+                                      * are nonzero. Data is organized as
+                                      * follows: if an entry is added to a
+                                      * row, it is first added to the @p cache
+                                      * variable, irrespective of whether an
+                                      * entry with same column number has
+                                      * already been added. Only if the cache
+                                      * is full do we flush it by removing
+                                      * duplicates, removing entries that are
+                                      * already stored in the @p entries
+                                      * array, sorting everything, and merging
+                                      * the two arrays.
+                                      *
+                                      * The reasoning behind this scheme is
+                                      * that memory allocation is expensive,
+                                      * and we only want to do it when really
+                                      * necessary. Previously (in deal.II
+                                      * versions up to 5.0), we used to store
+                                      * the column indices inside a std::set,
+                                      * but this would allocate 20 bytes each
+                                      * time we added an entry. Using the
+                                      * present scheme, we only need to
+                                      * allocate memory once for every 8 added
+                                      * entries, and we waste a lot less
+                                      * memory by not using a balanced tree
+                                      * for storing column indices.
+                                      *
+                                      * Since some functions that are @p const
+                                      * need to access the data of this
+                                      * object, but need to flush caches
+                                      * before, the @p flush_cache function is
+                                      * marked const, and the data members are
+                                      * marked @p mutable.
+                                      */
+    struct Line
+    {
+                                         /**
+                                          * Storage for the column indices of
+                                          * this row, unless they are still in
+                                          * the cache. This array is always
+                                          * kept sorted.
+                                          */
+        mutable std::vector<unsigned int> entries;
+
+                                         /**
+                                          * Size of the cache.
+                                          */
+        static const unsigned int cache_size = 8;
+        
+                                         /**
+                                          * Cache of entries that have not yet
+                                          * been written to @p entries;
+                                          */
+        mutable unsigned int cache[cache_size];
+
+                                         /**
+                                          * Number of entries in the cache.
+                                          */
+        mutable unsigned int cache_entries;
+
+                                         /**
+                                          * Constructor.
+                                          */
+        Line ();
+
+                                         /**
+                                          * Add the given column number to
+                                          * this line.
+                                          */
+        void add (const unsigned int col_num);
+
+                                         /**
+                                          * Flush the cache my merging it with
+                                          * the @p entries array.
+                                          */
+        void flush_cache () const;
+    };
+    
+        
 				     /**
 				      * Actual data: store for each
 				      * row the set of nonzero
 				      * entries.
 				      */
-    std::vector<std::set<unsigned int> > column_indices;
+    std::vector<Line> lines;
 };
 
 /*@}*/
@@ -367,6 +448,7 @@ CompressedSparsityPattern::n_rows () const
 }
 
 
+
 inline
 unsigned int
 CompressedSparsityPattern::n_cols () const
@@ -374,4 +456,226 @@ CompressedSparsityPattern::n_cols () const
   return cols;
 }
 
+
+
+inline
+void
+CompressedSparsityPattern::add (const unsigned int i,
+				const unsigned int j)
+{
+  Assert (i<rows, ExcInvalidIndex(i,rows));
+  Assert (j<cols, ExcInvalidIndex(j,cols));
+
+  lines[i].add (j);
+}
+
+
+
+inline
+CompressedSparsityPattern::Line::Line ()
+                :
+                cache_entries (0)
+{}
+
+
+
+inline
+unsigned int
+CompressedSparsityPattern::row_length (const unsigned int row) const
+{
+  Assert (row < n_rows(), ExcIndexRange (row, 0, n_rows()));
+  
+  lines[row].flush_cache ();
+  return lines[row].entries.size();
+}
+
+
+
+inline
+unsigned int
+CompressedSparsityPattern::column_number (const unsigned int row,
+					  const unsigned int index) const
+{
+  Assert (row < n_rows(), ExcIndexRange (row, 0, n_rows()));
+  Assert (index < lines[row].entries.size(),
+	  ExcIndexRange (index, 0, lines[row].entries.size()));
+
+  lines[row].flush_cache ();
+  return lines[row].entries[index];
+}
+
+
+
+inline
+void
+CompressedSparsityPattern::Line::add (const unsigned int j)
+{
+                                   // first check whether this entry is
+                                   // already in the cache. if so, we can
+                                   // safely return
+  for (unsigned int i=0; i<cache_entries; ++i)
+    if (cache[i] == j)
+      return;
+
+                                   // if not, see whether there is still some
+                                   // space in the cache. if not, then flush
+                                   // the cache first
+  if (cache_entries == cache_size)
+    flush_cache ();
+  
+  cache[cache_entries] = j;
+  ++cache_entries;
+}
+
+
+
+inline
+void
+CompressedSparsityPattern::Line::flush_cache () const
+{
+                                   // do nothing if the cache is empty
+  if (cache_entries == 0)
+    return;
+  
+                                   // first sort the entries in the cache, so
+                                   // that it is easier to merge it with the
+                                   // main array. note that due to the way
+                                   // add() inserts elements, there can be no
+                                   // duplicates in the cache
+                                   //
+                                   // do the sorting in a way that is fast for
+                                   // the small cache sizes we have
+                                   // here. basically, use bubble sort
+  switch (cache_entries)
+    {
+      case 1:
+      {
+        break;
+      }
+
+      case 2:
+      {
+        if (cache[1] < cache[0])
+          std::swap (cache[0], cache[1]);
+        break;
+      }
+
+      case 3:
+      {
+        if (cache[1] < cache[0])
+          std::swap (cache[0], cache[1]);
+        if (cache[2] < cache[1])
+          std::swap (cache[1], cache[2]);
+        if (cache[1] < cache[0])
+          std::swap (cache[0], cache[1]);
+        break;
+      }
+
+      case 4:
+      case 5:
+      case 6:
+      case 7:
+      {
+        for (unsigned int i=0; i<cache_entries; ++i)
+          for (unsigned int j=i+1; j<cache_entries; ++j)
+            if (cache[j] < cache[i])
+              std::swap (cache[i], cache[j]);
+        break;
+      }
+
+      default:
+      {
+        std::sort (&cache[0], &cache[cache_entries]);
+        break;
+      }
+    }
+
+                                   // next job is to merge the two
+                                   // arrays. special case the case that the
+                                   // original array is empty.
+  if (entries.size() == 0)
+    {
+      entries.resize (cache_entries);
+      for (unsigned int i=0; i<cache_entries; ++i)
+        entries[i] = cache[i];
+    }
+  else
+    {
+                                       // first count how many of the cache
+                                       // entries are already in the main
+                                       // array, so that we can efficiently
+                                       // allocate memory
+      unsigned int n_new_entries = 0;
+      {
+        unsigned int cache_position = 0;
+        unsigned int entry_position = 0;
+        while ((entry_position<entries.size()) &&
+               (cache_position<cache_entries))
+          {
+            ++n_new_entries;
+            if (entries[entry_position] < cache[cache_position])
+              ++entry_position;
+            else if (entries[entry_position] == cache[cache_position])
+              {
+                ++entry_position;
+                ++cache_position;
+              }
+            else
+              ++cache_position;
+          }
+
+                                         // scoop up leftovers in arrays
+        n_new_entries += (entries.size() - entry_position) +
+                         (cache_entries - cache_position);
+      }
+
+                                       // then allocate new memory and merge
+                                       // arrays
+      std::vector<unsigned int> new_entries;
+      new_entries.reserve (n_new_entries);
+      unsigned int cache_position = 0;
+      unsigned int entry_position = 0;
+      while ((entry_position<entries.size()) && (cache_position<cache_entries))
+        if (entries[entry_position] < cache[cache_position])
+          {
+            new_entries.push_back (entries[entry_position]);
+            ++entry_position;
+          }
+        else if (entries[entry_position] == cache[cache_position])
+          {
+            new_entries.push_back (entries[entry_position]);
+            ++entry_position;
+            ++cache_position;
+          }
+        else
+          {
+            new_entries.push_back (cache[cache_position]);
+            ++cache_position;
+          }
+
+                                       // copy remaining elements from the
+                                       // array that we haven't finished. note
+                                       // that at most one of the following
+                                       // loops will run at all
+      for (; entry_position < entries.size(); ++entry_position)
+        new_entries.push_back (entries[entry_position]);
+      for (; cache_position < cache_entries; ++cache_position)
+        new_entries.push_back (cache[cache_position]);
+      
+      Assert (new_entries.size() == n_new_entries,
+              ExcInternalError());
+      Assert (new_entries.capacity() == n_new_entries,
+              ExcInternalError());
+
+                                       // finally swap old and new array, and
+                                       // set cache size to zero
+      new_entries.swap (entries);
+    }
+  
+  cache_entries = 0;
+}
+
+
+
+
 #endif
diff --git a/deal.II/lac/source/compressed_sparsity_pattern.cc b/deal.II/lac/source/compressed_sparsity_pattern.cc
index 8536379115..806316172c 100644
--- a/deal.II/lac/source/compressed_sparsity_pattern.cc
+++ b/deal.II/lac/source/compressed_sparsity_pattern.cc
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2001, 2002, 2003 by the deal.II authors
+//    Copyright (C) 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
@@ -21,6 +21,11 @@
 #include <numeric>
 #include <functional>
 
+
+const unsigned int CompressedSparsityPattern::Line::cache_size;
+
+
+
 CompressedSparsityPattern::CompressedSparsityPattern ()
                 :
 		rows(0),
@@ -84,8 +89,8 @@ CompressedSparsityPattern::reinit (const unsigned int m,
   rows = m;
   cols = n;
 
-  std::vector<std::set<unsigned int> > new_column_indices (rows);
-  column_indices.swap (new_column_indices);
+  std::vector<Line> new_lines (rows);
+  lines.swap (new_lines);
 }
 
 
@@ -108,35 +113,28 @@ unsigned int
 CompressedSparsityPattern::max_entries_per_row () const
 {
   unsigned int m = 0;
-  for (unsigned int i=1; i<rows; ++i)
-    m = std::max (m, static_cast<unsigned int>(column_indices[i].size()));
-
+  for (unsigned int i=0; i<rows; ++i)
+    {
+      lines[i].flush_cache ();
+      m = std::max (m, static_cast<unsigned int>(lines[i].entries.size()));
+    }
+  
   return m;
 }
 
 
 
-void
-CompressedSparsityPattern::add (const unsigned int i,
-				const unsigned int j)
-{
-  Assert (i<rows, ExcInvalidIndex(i,rows));
-  Assert (j<cols, ExcInvalidIndex(j,cols));
-
-				   // the std::set automatically
-				   // assures uniqueness and
-				   // sortedness of the column indices
-  column_indices[i].insert (j);
-}
-
-
 bool 
 CompressedSparsityPattern::exists (const unsigned int i,
 				   const unsigned int j) const
 {
   Assert (i<rows, ExcInvalidIndex(i,rows));
   Assert (j<cols, ExcInvalidIndex(j,cols));
-  return (column_indices[i].find(j) != column_indices[i].end());
+
+  lines[i].flush_cache();
+  return std::binary_search (lines[i].entries.begin(),
+                             lines[i].entries.end(),
+                             j);
 }
 
 
@@ -158,67 +156,59 @@ CompressedSparsityPattern::symmetrize ()
 				   // be called on elements that
 				   // already exist without any harm
   for (unsigned int row=0; row<rows; ++row)
-    for (std::set<unsigned int>::const_iterator i=column_indices[row].begin();
-	 i!=column_indices[row].end(); ++i)
+    {
+      lines[row].flush_cache ();
+      for (std::vector<unsigned int>::const_iterator
+             j=lines[row].entries.begin();
+           j != lines[row].entries.end();
+           ++j)
 				       // add the transpose entry if
 				       // this is not the diagonal
-      if (row != *i)
-	add (*i, row);
+        if (row != *j)
+          add (*j, row);
+    }
 }
 
 
 
 void
 CompressedSparsityPattern::print_gnuplot (std::ostream &out) const
-{
+{ 
   for (unsigned int row=0; row<rows; ++row)
-    for (std::set<unsigned int>::const_iterator i=column_indices[row].begin();
-	 i!=column_indices[row].end(); ++i)
-				       // while matrix entries are
-				       // usually written (i,j),
-				       // with i vertical and j
-				       // horizontal, gnuplot output
-				       // is x-y, that is we have to
-				       // exchange the order of
-				       // output
-      out << *i << " " << -static_cast<signed int>(row) << std::endl;
+    {
+      lines[row].flush_cache ();
+      for (std::vector<unsigned int>::const_iterator
+             j=lines[row].entries.begin();
+           j != lines[row].entries.end(); ++j)
+                                         // while matrix entries are usually
+                                         // written (i,j), with i vertical and
+                                         // j horizontal, gnuplot output is
+                                         // x-y, that is we have to exchange
+                                         // the order of output
+        out << *j << " " << -static_cast<signed int>(row) << std::endl;
+    }
+      
 
   AssertThrow (out, ExcIO());
 }
 
 
 
-unsigned int
-CompressedSparsityPattern::row_length (const unsigned int row) const
-{
-  return column_indices[row].size();
-}
-
-
-
-unsigned int
-CompressedSparsityPattern::column_number (const unsigned int row,
-					  const unsigned int index) const
-{
-  Assert (index < column_indices[row].size(),
-	  ExcIndexRange (index, 0, column_indices[row].size()));
-  std::set<unsigned int>::const_iterator p = column_indices[row].begin();
-  std::advance (p, index);
-  return *p;
-}
-
-
-
 unsigned int
 CompressedSparsityPattern::bandwidth () const
 {
   unsigned int b=0;
   for (unsigned int row=0; row<rows; ++row)
-    for (std::set<unsigned int>::const_iterator i=column_indices[row].begin();
-	 i!=column_indices[row].end(); ++i)
-      if (static_cast<unsigned int>(std::abs(static_cast<int>(row-*i))) > b)
-	b = std::abs(static_cast<signed int>(row-*i));
-
+    {
+      lines[row].flush_cache ();
+
+      for (std::vector<unsigned int>::const_iterator
+             j=lines[row].entries.begin();
+           j != lines[row].entries.end(); ++j)
+        if (static_cast<unsigned int>(std::abs(static_cast<int>(row-*j))) > b)
+          b = std::abs(static_cast<signed int>(row-*j));
+    }
+  
   return b;
 }
 
@@ -229,6 +219,10 @@ CompressedSparsityPattern::n_nonzero_elements () const
 {
   unsigned int n=0;
   for (unsigned int i=0; i<rows; ++i)
-    n += column_indices[i].size();
+    {
+      lines[i].flush_cache ();
+      n += lines[i].entries.size();
+    }
+  
   return n;
 }
-- 
2.39.5