From a222ab7c5c5f3d192ac23fd31bb2a20fa375de78 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 25 Apr 2001 16:10:02 +0000 Subject: [PATCH] Use a much faster replacement for std::lower_bound which took about one third of the time to perform the SSOR preconditioner (which in turn took about half the total time needed to run my wave program). git-svn-id: https://svn.dealii.org/trunk@4488 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/2001/c-3-1.html | 15 +++ deal.II/lac/include/lac/sparse_matrix.h | 2 +- .../lac/include/lac/sparse_matrix.templates.h | 16 +-- deal.II/lac/include/lac/sparsity_pattern.h | 115 ++++++++++++++++++ deal.II/lac/source/sparsity_pattern.cc | 29 ++--- 5 files changed, 155 insertions(+), 22 deletions(-) diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index ff3b6834cb..30f2e8e513 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -234,6 +234,21 @@ documentation, etc.

lac

    +
  1. + New: There is now a (private) function SparsityPattern::optimized_lower_bound + that is used as an optimized replacement for std::lower_bound 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. +
    + (WB 2001/04/25) +

    +
  2. New: There is now a function Vector::scale(Vector) diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index 52ddead5d4..c94aafe27c 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -876,7 +876,6 @@ class SparseMatrix : public Subscriptor const std::pair interval, somenumber *partial_norm) const; - // make all other sparse matrices // friends template friend class SparseMatrix; @@ -886,6 +885,7 @@ class SparseMatrix : public Subscriptor /*---------------------- Inline functions -----------------------------------*/ + template inline unsigned int SparseMatrix::m () const diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index 4d116df7ad..66b17ce955 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -32,6 +32,8 @@ #endif + + template SparseMatrix::SparseMatrix () : cols(0), @@ -774,12 +776,12 @@ SparseMatrix::precondition_SSOR (Vector& 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; jcolnums[j]); @@ -798,9 +800,9 @@ SparseMatrix::precondition_SSOR (Vector& 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(row)) - + = (SparsityPattern::optimized_lower_bound (&cols->colnums[*rowstart_ptr+1], + &cols->colnums[*(rowstart_ptr+1)], + static_cast(row)) - &cols->colnums[0]); for (unsigned int j=first_right_of_diagonal_index; j<*(rowstart_ptr+1); ++j) if (cols->colnums[j] > static_cast(row)) diff --git a/deal.II/lac/include/lac/sparsity_pattern.h b/deal.II/lac/include/lac/sparsity_pattern.h index 34caa5f12f..898dbe6fee 100644 --- a/deal.II/lac/include/lac/sparsity_pattern.h +++ b/deal.II/lac/include/lac/sparsity_pattern.h @@ -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 diff --git a/deal.II/lac/source/sparsity_pattern.cc b/deal.II/lac/source/sparsity_pattern.cc index dac608db76..8ab04cf613 100644 --- a/deal.II/lac/source/sparsity_pattern.cc +++ b/deal.II/lac/source/sparsity_pattern.cc @@ -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