*/
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;
};
/*@}*/
}
+
inline
unsigned int
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
// $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
#include <numeric>
#include <functional>
+
+const unsigned int CompressedSparsityPattern::Line::cache_size;
+
+
+
CompressedSparsityPattern::CompressedSparsityPattern ()
:
rows(0),
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);
}
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);
}
// 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;
}
{
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;
}
projects / dealii-svn.git / commitdiff