--- /dev/null
+//---------------------------------------------------------------------------
+// $Id: compressed_sparsity_pattern.h 14038 2006-10-23 02:46:34Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__compressed_set_sparsity_pattern_h
+#define __deal2__compressed_set_sparsity_pattern_h
+
+
+#include <base/config.h>
+#include <base/subscriptor.h>
+#include <lac/exceptions.h>
+
+#include <vector>
+#include <algorithm>
+#include <set>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <typename number> class SparseMatrix;
+
+/*! @addtogroup Matrix1
+ *@{
+ */
+
+
+/**
+ * This class acts as an intermediate form of the
+ * SparsityPattern class. From the interface it mostly
+ * represents a SparsityPattern object that is kept compressed
+ * at all times. However, since the final sparsity pattern is not
+ * known while constructing it, keeping the pattern compressed at all
+ * times can only be achieved at the expense of either increased
+ * memory or run time consumption upon use. The main purpose of this
+ * class is to avoid some memory bottlenecks, so we chose to implement
+ * it memory conservative, but the chosen data format is too unsuited
+ * to be used for actual matrices. It is therefore necessary to first
+ * copy the data of this object over to an object of type
+ * SparsityPattern before using it in actual matrices.
+ *
+ * Another viewpoint is that this class does not need up front
+ * allocation of a certain amount of memory, but grows as necessary.
+ *
+ *
+ * <h3>Rationale</h3>
+ *
+ * When constructing the sparsity pattern of a matrix, you usually
+ * first have to provide an empty sparsity pattern object with a fixed
+ * maximal number of entries per row. To find out about this maximal
+ * row length, one usually calls the function
+ * DoFHandler::max_couplings_between_dofs() which returns an estimate for
+ * that quantity or DoFTools::compute_row_length_vector(), which gives
+ * a better estimate for each row. While this estimate is usually
+ * quite good in 2d and exact in 1d, it is often significantly too
+ * large in 3d and especially for higher order elements. Furthermore,
+ * normally only a small fraction of the rows of a matrix will end up
+ * having the maximal number of nonzero entries per row (usually those
+ * nodes adjacent to hanging nodes), most have much less. In effect,
+ * the empty SparsityPattern object has allocated much too much
+ * memory. Although this unnecessarily allocated memory is later freed
+ * when SparsityPattern::compress() is called, this overallocation
+ * has, with higher order elements and in 3d, sometimes been so large
+ * that the program aborted due to lack of memory.
+ *
+ * This class therefore provides an alternative representation of a
+ * sparsity pattern: we don't specify a maximal row length initially,
+ * but store a set of column indices indicating possible nonzero
+ * entries in the sparsity pattern for each row. This is very much
+ * like the final "compressed" format used in the
+ * SparsityPattern object after compression, but uses a less
+ * compact memory storage format, since the exact number of entries
+ * per row is only known a posteriori and since it may change (for the
+ * SparsityPattern class, no more changes are allowed after
+ * compressing it). We can therefore not store all the column indices
+ * in a big array, but have to use a vector of sets. This can later be
+ * used to actually initialize a SparsityPattern object with the
+ * then final set of necessary indices.
+ *
+ *
+ * <h3>Interface</h3>
+ *
+ * Since this class is intended as an intermediate replacement of the
+ * SparsityPattern class, it has mostly the same interface, with
+ * small changes where necessary. In particular, the add()
+ * function, and the functions inquiring properties of the sparsity
+ * pattern are the same.
+ *
+ * <h3>Notes</h3>
+ *
+ * This class is a variation of the CompressedSparsityPattern class.
+ * Instead of using sorted vectors together with a caching algorithm
+ * for storing the column indices of NZ entries, the std::set
+ * container is used. This solution might not be the fastest in
+ * all situations, but seems to work much better than the
+ * CompressedSparsityPattern in the context of hp-adaptivity.
+ * On the other hand, a benchmark where NZ entries were randomly inserted
+ * into the sparsity pattern revealed that this class is slower
+ * by a factor 4-6 in this situation. Hence, currently the suggestion
+ * is to carefully analyse which of the CompressedSparsityPattern
+ * classes works best in a certain setting. An algorithm which
+ * performs equally well in all situations still has to be found.
+ *
+ * <h3>Usage</h3>
+ *
+ * Use this class as follows:
+ * @verbatim
+ * CompressedSetSparsityPattern compressed_pattern (dof_handler.n_dofs());
+ * DoFTools::make_sparsity_pattern (dof_handler,
+ * compressed_pattern);
+ * constraints.condense (compressed_pattern);
+ *
+ * SparsityPattern sp;
+ * sp.copy_from (compressed_pattern);
+ * @endverbatim
+ *
+ *
+ * @author Wolfgang Bangerth, 2001
+ */
+class CompressedSetSparsityPattern : public Subscriptor
+{
+ public:
+ typedef std::set<unsigned int>::iterator CSSPIterator;
+
+
+ /**
+ * Initialize the matrix empty,
+ * that is with no memory
+ * allocated. This is useful if
+ * you want such objects as
+ * member variables in other
+ * classes. You can make the
+ * structure usable by calling
+ * the reinit() function.
+ */
+ CompressedSetSparsityPattern ();
+
+ /**
+ * Copy constructor. This constructor is
+ * only allowed to be called if the
+ * matrix structure to be copied is
+ * empty. This is so in order to prevent
+ * involuntary copies of objects for
+ * temporaries, which can use large
+ * amounts of computing time. However,
+ * copy constructors are needed if yo
+ * want to use the STL data types on
+ * classes like this, e.g. to write such
+ * statements like <tt>v.push_back
+ * (CompressedSetSparsityPattern());</tt>,
+ * with @p v a vector of @p
+ * CompressedSetSparsityPattern objects.
+ *
+ * Usually, it is sufficient to
+ * use the explicit keyword to
+ * disallow unwanted temporaries,
+ * but for the STL vectors, this
+ * does not work. Since copying a
+ * structure like this is not
+ * useful anyway because multiple
+ * matrices can use the same
+ * sparsity structure, copies are
+ * only allowed for empty
+ * objects, as described above.
+ */
+ CompressedSetSparsityPattern (const CompressedSetSparsityPattern &);
+
+ /**
+ * Initialize a rectangular
+ * matrix with @p m rows and
+ * @p n columns.
+ */
+ CompressedSetSparsityPattern (const unsigned int m,
+ const unsigned int n);
+
+ /**
+ * Initialize a square matrix of
+ * dimension @p n.
+ */
+ CompressedSetSparsityPattern (const unsigned int n);
+
+ /**
+ * Copy operator. For this the
+ * same holds as for the copy
+ * constructor: it is declared,
+ * defined and fine to be called,
+ * but the latter only for empty
+ * objects.
+ */
+ CompressedSetSparsityPattern & operator = (const CompressedSetSparsityPattern &);
+
+ /**
+ * Reallocate memory and set up
+ * data structures for a new
+ * matrix with @p m rows and
+ * @p n columns, with at most
+ * max_entries_per_row() nonzero
+ * entries per row.
+ */
+ void reinit (const unsigned int m,
+ const unsigned int n);
+
+ /**
+ * Since this object is kept
+ * compressed at all times anway,
+ * this function does nothing,
+ * but is declared to make the
+ * interface of this class as
+ * much alike as that of the
+ * SparsityPattern class.
+ */
+ void compress ();
+
+ /**
+ * Return whether the object is
+ * empty. It is empty if no
+ * memory is allocated, which is
+ * the same as that both
+ * dimensions are zero.
+ */
+ bool empty () const;
+
+ /**
+ * Return the maximum number of
+ * entries per row. Note that
+ * this number may change as
+ * entries are added.
+ */
+ unsigned int max_entries_per_row () const;
+
+ /**
+ * Add a nonzero entry to the
+ * matrix. If the entry already
+ * exists, nothing bad happens.
+ */
+ void add (const unsigned int i,
+ const unsigned int j);
+
+ /**
+ * Check if a value at a certain
+ * position may be non-zero.
+ */
+ bool exists (const unsigned int i,
+ const unsigned int j) const;
+
+ /**
+ * Make the sparsity pattern
+ * symmetric by adding the
+ * sparsity pattern of the
+ * transpose object.
+ *
+ * This function throws an
+ * exception if the sparsity
+ * pattern does not represent a
+ * square matrix.
+ */
+ void symmetrize ();
+
+ /**
+ * Print the sparsity of the matrix in a
+ * format that @p gnuplot understands and
+ * which can be used to plot the sparsity
+ * pattern in a graphical way. The format
+ * consists of pairs <tt>i j</tt> of
+ * nonzero elements, each representing
+ * one entry of this matrix, one per line
+ * of the output file. Indices are
+ * counted from zero on, as usual. Since
+ * sparsity patterns are printed in the
+ * same way as matrices are displayed, we
+ * print the negative of the column
+ * index, which means that the
+ * <tt>(0,0)</tt> element is in the top
+ * left rather than in the bottom left
+ * corner.
+ *
+ * Print the sparsity pattern in
+ * gnuplot by setting the data style
+ * to dots or points and use the
+ * @p plot command.
+ */
+ void print_gnuplot (std::ostream &out) const;
+
+ /**
+ * Return number of rows of this
+ * matrix, which equals the dimension
+ * of the image space.
+ */
+ unsigned int n_rows () const;
+
+ /**
+ * Return number of columns of this
+ * matrix, which equals the dimension
+ * of the range space.
+ */
+ unsigned int n_cols () const;
+
+ /**
+ * Number of entries in a specific row.
+ */
+ unsigned int row_length (const unsigned int row) const;
+
+ /**
+ * Access to column number field.
+ * Return the column number of
+ * the @p indexth entry in @p row.
+ */
+ CSSPIterator row_begin (const unsigned int row) const;
+ CSSPIterator row_end (const unsigned int row) const;
+
+
+ /**
+ * Compute the bandwidth of the matrix
+ * represented by this structure. The
+ * bandwidth is the maximum of
+ * $|i-j|$ for which the index pair
+ * $(i,j)$ represents a nonzero entry
+ * of the matrix.
+ */
+ unsigned int bandwidth () const;
+
+ /**
+ * Return the number of nonzero elements
+ * of this matrix. Actually, it returns
+ * the number of entries in the sparsity
+ * pattern; if any of the entries should
+ * happen to be zero, it is counted
+ * anyway.
+ *
+ * This function may only be called if
+ * the matrix struct is compressed. It
+ * does not make too much sense otherwise
+ * anyway.
+ */
+ unsigned int n_nonzero_elements () const;
+
+ private:
+ /**
+ * Number of rows that this sparsity
+ * structure shall represent.
+ */
+ unsigned int rows;
+
+ /**
+ * Number of columns that this sparsity
+ * structure shall represent.
+ */
+ 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 #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 flush_cache() function is
+ * marked const, and the data members are
+ * marked @p mutable.
+ *
+ * A small testseries about the size of
+ * the cache showed that the run time of
+ * a small program just testing the
+ * compressed sparsity pattern element
+ * insertion routine ran for 3.6 seconds
+ * with a cache size of 8, and 4.2
+ * seconds with a cache size of 16. We
+ * deem even smaller cache sizes
+ * undesirable, since they lead to more
+ * memory allocations, while larger cache
+ * sizes lead to waste of memory. The
+ * original version of this class, with
+ * one std::set per row took 8.2 seconds
+ * on the same program.
+ */
+ struct Line
+ {
+ std::set <unsigned int> entries;
+
+ /**
+ * Constructor.
+ */
+ Line ();
+
+ /**
+ * Add the given column number to
+ * this line.
+ */
+ void add (const unsigned int col_num);
+ };
+
+
+ /**
+ * Actual data: store for each
+ * row the set of nonzero
+ * entries.
+ */
+ std::vector<Line> lines;
+};
+
+/*@}*/
+/*---------------------- Inline functions -----------------------------------*/
+
+
+inline
+void
+CompressedSetSparsityPattern::Line::add (const unsigned int j)
+{
+ entries.insert (j);
+}
+
+
+
+inline
+unsigned int
+CompressedSetSparsityPattern::n_rows () const
+{
+ return rows;
+}
+
+
+
+inline
+unsigned int
+CompressedSetSparsityPattern::n_cols () const
+{
+ return cols;
+}
+
+
+
+inline
+void
+CompressedSetSparsityPattern::add (const unsigned int i,
+ const unsigned int j)
+{
+ Assert (i<rows, ExcIndexRange(i, 0, rows));
+ Assert (j<cols, ExcIndexRange(j, 0, cols));
+
+ lines[i].add (j);
+}
+
+
+
+inline
+CompressedSetSparsityPattern::Line::Line ()
+{}
+
+
+
+inline
+unsigned int
+CompressedSetSparsityPattern::row_length (const unsigned int row) const
+{
+ Assert (row < n_rows(), ExcIndexRange (row, 0, n_rows()));
+
+ return lines[row].entries.size();
+}
+
+
+inline
+CompressedSetSparsityPattern::CSSPIterator
+CompressedSetSparsityPattern::row_begin (const unsigned int row) const
+{
+ return (lines[row].entries.begin ());
+}
+
+
+inline
+CompressedSetSparsityPattern::CSSPIterator
+CompressedSetSparsityPattern::row_end (const unsigned int row) const
+{
+ return (lines[row].entries.end ());
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id: compressed_sparsity_pattern.cc 14038 2006-10-23 02:46:34Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+#include <lac/compressed_set_sparsity_pattern.h>
+
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+#include <functional>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+CompressedSetSparsityPattern::CompressedSetSparsityPattern ()
+ :
+ rows(0),
+ cols(0)
+{}
+
+
+
+CompressedSetSparsityPattern::
+CompressedSetSparsityPattern (const CompressedSetSparsityPattern &s)
+ :
+ Subscriptor(),
+ rows(0),
+ cols(0)
+{
+ Assert (s.rows == 0, ExcInvalidConstructorCall());
+ Assert (s.cols == 0, ExcInvalidConstructorCall());
+}
+
+
+
+CompressedSetSparsityPattern::CompressedSetSparsityPattern (const unsigned int m,
+ const unsigned int n)
+ :
+ rows(0),
+ cols(0)
+{
+ reinit (m,n);
+}
+
+
+
+CompressedSetSparsityPattern::CompressedSetSparsityPattern (const unsigned int n)
+ :
+ rows(0),
+ cols(0)
+{
+ reinit (n,n);
+}
+
+
+
+CompressedSetSparsityPattern &
+CompressedSetSparsityPattern::operator = (const CompressedSetSparsityPattern &s)
+{
+ Assert (s.rows == 0, ExcInvalidConstructorCall());
+ Assert (s.cols == 0, ExcInvalidConstructorCall());
+
+ Assert (rows == 0, ExcInvalidConstructorCall());
+ Assert (cols == 0, ExcInvalidConstructorCall());
+
+ return *this;
+}
+
+
+
+void
+CompressedSetSparsityPattern::reinit (const unsigned int m,
+ const unsigned int n)
+{
+ rows = m;
+ cols = n;
+
+ std::vector<Line> new_lines (rows);
+ lines.swap (new_lines);
+}
+
+
+
+void
+CompressedSetSparsityPattern::compress ()
+{}
+
+
+
+bool
+CompressedSetSparsityPattern::empty () const
+{
+ return ((rows==0) && (cols==0));
+}
+
+
+
+unsigned int
+CompressedSetSparsityPattern::max_entries_per_row () const
+{
+ unsigned int m = 0;
+ for (unsigned int i=0; i<rows; ++i)
+ {
+ m = std::max (m, static_cast<unsigned int>(lines[i].entries.size()));
+ }
+
+ return m;
+}
+
+
+
+bool
+CompressedSetSparsityPattern::exists (const unsigned int i,
+ const unsigned int j) const
+{
+ Assert (i<rows, ExcIndexRange(i, 0, rows));
+ Assert (j<cols, ExcIndexRange(j, 0, cols));
+
+ return (lines[i].entries.find (j) != lines[i].entries.end ());
+}
+
+
+
+void
+CompressedSetSparsityPattern::symmetrize ()
+{
+ Assert (rows==cols, ExcNotQuadratic());
+
+ // loop over all elements presently
+ // in the sparsity pattern and add
+ // the transpose element. note:
+ //
+ // 1. that the sparsity pattern
+ // changes which we work on, but
+ // not the present row
+ //
+ // 2. that the @p{add} function can
+ // 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
+ j=lines[row].entries.begin();
+ j != lines[row].entries.end();
+ ++j)
+ // add the transpose entry if
+ // this is not the diagonal
+ if (row != *j)
+ add (*j, row);
+ }
+}
+
+
+
+void
+CompressedSetSparsityPattern::print_gnuplot (std::ostream &out) const
+{
+ for (unsigned int row=0; row<rows; ++row)
+ {
+ for (std::set<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
+CompressedSetSparsityPattern::bandwidth () const
+{
+ unsigned int b=0;
+ for (unsigned int row=0; row<rows; ++row)
+ {
+ for (std::set<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
+CompressedSetSparsityPattern::n_nonzero_elements () const
+{
+ unsigned int n=0;
+ for (unsigned int i=0; i<rows; ++i)
+ {
+ n += lines[i].entries.size();
+ }
+
+ return n;
+}
+
+DEAL_II_NAMESPACE_CLOSE
projects / dealii-svn.git / commitdiff