]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Use a much faster replacement for std::lower_bound which took about one third of...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 25 Apr 2001 16:10:02 +0000 (16:10 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 25 Apr 2001 16:10:02 +0000 (16:10 +0000)
git-svn-id: https://svn.dealii.org/trunk@4488 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/2001/c-3-1.html
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix.templates.h
deal.II/lac/include/lac/sparsity_pattern.h
deal.II/lac/source/sparsity_pattern.cc

index ff3b6834cbd653dd2974680a55cbe4d80e1873a7..30f2e8e513303acfe2b0e77494048b822090f91d 100644 (file)
@@ -234,6 +234,21 @@ documentation, etc</a>.
 <h3>lac</h3>
 
 <ol>
+  <li> <p>
+       New: There is now a (private) function <code
+       class="member">SparsityPattern::optimized_lower_bound</code>
+       that is used as an optimized replacement for <code
+       class="member">std::lower_bound</code> for searching in the
+       column number arrays. It unrolls small loops and it also seems
+       that the compiler is able to optimized it better due to
+       eliminated template parameters, making it about twice as fast
+       as the standard implementation. In effect, it also speeds up
+       the SSOR preconditioner that spent about one third of its time
+       in that function by approximately 15 per cent.
+       <br>
+       (WB 2001/04/25)
+       </p>
+
   <li> <p>
        New: There is now a function <code
        class="member">Vector::scale(Vector)</code>
index 52ddead5d46c9d2e13ffa8056212603701a8232d..c94aafe27c62cac889dd4d98173159aabc62513f 100644 (file)
@@ -876,7 +876,6 @@ class SparseMatrix : public Subscriptor
                            const std::pair<unsigned int,unsigned int> interval,
                            somenumber               *partial_norm) const;
 
-
                                     // make all other sparse matrices
                                     // friends
     template <typename somenumber> friend class SparseMatrix;
@@ -886,6 +885,7 @@ class SparseMatrix : public Subscriptor
 /*---------------------- Inline functions -----------------------------------*/
 
 
+
 template <typename number>
 inline
 unsigned int SparseMatrix<number>::m () const
index 4d116df7ad440633db65dc9db6b04ab7f9b26b7a..66b17ce955bdd119bc7f60353c53f38025e1ff9d 100644 (file)
@@ -32,6 +32,8 @@
 #endif
 
 
+
+
 template <typename number>
 SparseMatrix<number>::SparseMatrix () :
                cols(0),
@@ -774,12 +776,12 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>& dst,
                                       // line denotes the diagonal element,
                                       // which we need not check.
       const unsigned int first_right_of_diagonal_index
-       = (std::lower_bound (&cols->colnums[*rowstart_ptr+1],
-                            &cols->colnums[*(rowstart_ptr+1)],
-                            row)
+       = (SparsityPattern::optimized_lower_bound (&cols->colnums[*rowstart_ptr+1],
+                                                  &cols->colnums[*(rowstart_ptr+1)],
+                                                  row)
           -
           &cols->colnums[0]);
-                                      
+      
       for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
        *dst_ptr -= om* val[j] * dst(cols->colnums[j]);
 
@@ -798,9 +800,9 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>& dst,
   for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr)
     {
       const unsigned int first_right_of_diagonal_index
-       = (std::lower_bound (&cols->colnums[*rowstart_ptr+1],
-                            &cols->colnums[*(rowstart_ptr+1)],
-                            static_cast<unsigned int>(row)) -
+       = (SparsityPattern::optimized_lower_bound (&cols->colnums[*rowstart_ptr+1],
+                                                  &cols->colnums[*(rowstart_ptr+1)],
+                                                  static_cast<unsigned int>(row)) -
           &cols->colnums[0]);
       for (unsigned int j=first_right_of_diagonal_index; j<*(rowstart_ptr+1); ++j)
        if (cols->colnums[j] > static_cast<unsigned int>(row))
index 34caa5f12ffa452154c8e42c7ee7b4669f626ba1..898dbe6feeedeb57d86fe4c732f41e6b0f083788 100644 (file)
@@ -690,6 +690,28 @@ class SparsityPattern : public Subscriptor
                                      */
     bool compressed;
 
+                                    /**
+                                     * Optimized replacement for
+                                     * @p{std::lower_bound} for
+                                     * searching within the range of
+                                     * column indices. Slashes
+                                     * execution time by
+                                     * approximately one half for the
+                                     * present application, partly
+                                     * because we have eliminated
+                                     * templates and the compiler
+                                     * seems to be able to optimize
+                                     * better, and partly because the
+                                     * binary search is replaced by a
+                                     * linear search for small loop
+                                     * lengths.
+                                     */
+    static
+    const unsigned int * const
+    optimized_lower_bound (const unsigned int *first,
+                          const unsigned int *last,
+                          const unsigned int &val);
+    
                                     /**
                                      * Make all sparse matrices
                                      * friends of this class.
@@ -701,6 +723,99 @@ class SparsityPattern : public Subscriptor
 /*---------------------- Inline functions -----------------------------------*/
 
 
+inline
+const unsigned int * const
+SparsityPattern::optimized_lower_bound (const unsigned int *first,
+                                       const unsigned int *last,
+                                       const unsigned int &val)
+{
+                                  // this function is mostly copied
+                                  // over from the STL __lower_bound
+                                  // function, but with template args
+                                  // replaced by the actual data
+                                  // types needed here, and above all
+                                  // with a rolled out search on the
+                                  // last four elements
+  unsigned int len = last-first;
+
+  if (len==0)
+    return first;
+  
+  while (true)
+    {
+                                      // if length equals 8 or less,
+                                      // then do a rolled out
+                                      // search. use a switch without
+                                      // breaks for that and roll-out
+                                      // the loop somehow
+      if (len < 8)
+       {
+         switch (len)
+           {
+             case 7:
+                   if (*first >= val)
+                     return first;
+                   ++first;
+             case 6:
+                   if (*first >= val)
+                     return first;
+                   ++first;
+             case 5:
+                   if (*first >= val)
+                     return first;
+                   ++first;
+             case 4:
+                   if (*first >= val)
+                     return first;
+                   ++first;
+             case 3:
+                   if (*first >= val)
+                     return first;
+                   ++first;
+             case 2:
+                   if (*first >= val)
+                     return first;
+                   ++first;
+             case 1:
+                   if (*first >= val)
+                     return first;
+                   return first+1;
+             default:
+                                                    // indices seem
+                                                    // to not be
+                                                    // sorted
+                                                    // correctly!? or
+                                                    // did len
+                                                    // become==0
+                                                    // somehow? that
+                                                    // shouln't have
+                                                    // happened
+                   Assert (false, ExcInternalError());
+           };
+       };
+      
+      
+      
+      const unsigned int         half   = len >> 1;
+      const unsigned int * const middle = first + half;
+      
+                                      // if the value is larger than
+                                      // that pointed to by the
+                                      // middle pointer, then the
+                                      // insertion point must be
+                                      // right of it
+      if (*middle < val)
+       {
+         first = middle + 1;
+         len  -= half + 1;
+       }
+      else
+       len = half;
+    }
+}
+
+
+
 inline
 unsigned int
 SparsityPattern::n_rows () const
index dac608db76b6deb020c86e63e23473ee0e33fc19..8ab04cf613999d264779ac0bcd28029255d383a3 100644 (file)
@@ -496,6 +496,7 @@ SparsityPattern::max_entries_per_row () const
 
 
 
+inline
 unsigned int
 SparsityPattern::operator () (const unsigned int i,
                              const unsigned int j) const
@@ -505,19 +506,16 @@ SparsityPattern::operator () (const unsigned int i,
   Assert (j<cols, ExcInvalidIndex(j,cols));
   Assert (compressed, ExcNotCompressed());
 
+                                  // let's see whether there is
+                                  // something in this line
+  if (rowstart[i] == rowstart[i+1])
+    return invalid_entry;
+  
                                   // check first entry separately, since
                                   // for square matrices this is
-                                  // the diagonal entry (check only
-                                  // if a first entry exists)
-  if (rowstart[i] != rowstart[i+1]) 
-    {
-      if (j == colnums[rowstart[i]])
-       return rowstart[i];
-    }
-  else
-                                    // no first entry exists for this
-                                    // line
-    return invalid_entry;
+                                  // the diagonal entry
+  if ((i==j) && (rows==cols))
+    return rowstart[i];
 
                                   // all other entries are sorted, so
                                   // we can use a binary seach algorithm
@@ -529,9 +527,12 @@ SparsityPattern::operator () (const unsigned int i,
                                   // at the top of this function, so it
                                   // may not be called for noncompressed
                                   // structures.
-  const unsigned int * const p = std::lower_bound (&colnums[rowstart[i]+1],
-                                                  &colnums[rowstart[i+1]],
-                                                  j);
+  const unsigned int * const sorted_region_start = (rows==cols ?
+                                                   &colnums[rowstart[i]+1] :
+                                                   &colnums[rowstart[i]]);
+  const unsigned int * const p = optimized_lower_bound (sorted_region_start,
+                                                       &colnums[rowstart[i+1]],
+                                                       j);
   if ((*p == j) &&
       (p != &colnums[rowstart[i+1]]))
     return (p - &colnums[0]);

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.