]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added an alternative implementation of a compressed
authorkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Jun 2007 15:22:40 +0000 (15:22 +0000)
committerkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Jun 2007 15:22:40 +0000 (15:22 +0000)
sparsity pattern, which is based on the std::set class.
Due to unknown reasons, this class seems to perform
better than the original CSP in the context of
hp-adaptivity.

git-svn-id: https://svn.dealii.org/trunk@14758 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/compressed_set_sparsity_pattern.h [new file with mode: 0644]
deal.II/lac/include/lac/sparsity_pattern.h
deal.II/lac/source/compressed_set_sparsity_pattern.cc [new file with mode: 0644]
deal.II/lac/source/sparsity_pattern.cc

diff --git a/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h
new file mode 100644 (file)
index 0000000..63cb2fb
--- /dev/null
@@ -0,0 +1,510 @@
+//---------------------------------------------------------------------------
+//    $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
index 0c15457dcfdc5c393c6695057c8479d1fcbad0a1..9af439f4d094b2d3d8e7601d3f815eb9704b3e90 100644 (file)
@@ -29,6 +29,7 @@ template <typename number> class SparseMatrix;
 template <class VECTOR> class VectorSlice;
 
 class CompressedSparsityPattern;
+class CompressedSetSparsityPattern;
 
 
 
@@ -968,6 +969,20 @@ class SparsityPattern : public Subscriptor
     void copy_from (const CompressedSparsityPattern &csp,
                    const bool optimize_diagonal = true);
 
+
+                                    /**
+                                     * Copy data from an object of
+                                     * type
+                                     * CompressedSetSparsityPattern.
+                                     * Previous content of this
+                                     * object is lost, and the
+                                     * sparsity pattern is in
+                                     * compressed mode afterwards.
+                                     */
+    void copy_from (const CompressedSetSparsityPattern &csp,
+                   const bool optimize_diagonal = true);
+
+
                                     /**
                                      * Take a full matrix and use its
                                      * nonzero entries to generate a
diff --git a/deal.II/lac/source/compressed_set_sparsity_pattern.cc b/deal.II/lac/source/compressed_set_sparsity_pattern.cc
new file mode 100644 (file)
index 0000000..fd50f31
--- /dev/null
@@ -0,0 +1,220 @@
+//---------------------------------------------------------------------------
+//    $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
index ca058b9deaee728a73f80eebbc13c52709122bc7..05e1d4f3481b87d3418c215ea963b6d389f3cb07 100644 (file)
@@ -15,6 +15,7 @@
 #include <lac/sparsity_pattern.h>
 #include <lac/full_matrix.h>
 #include <lac/compressed_sparsity_pattern.h>
+#include <lac/compressed_set_sparsity_pattern.h>
 
 #include <iostream>
 #include <iomanip>
@@ -597,6 +598,57 @@ SparsityPattern::copy_from (const CompressedSparsityPattern &csp,
 }
 
 
+
+void
+SparsityPattern::copy_from (const CompressedSetSparsityPattern &csp,
+                           const bool optimize_diag) 
+{
+                                  // first determine row lengths for
+                                  // each row. if the matrix is
+                                  // quadratic, then we might have to
+                                  // add an additional entry for the
+                                  // diagonal, if that is not yet
+                                  // present. as we have to call
+                                  // compress anyway later on, don't
+                                  // bother to check whether that
+                                  // diagonal entry is in a certain
+                                  // row or not
+  const bool is_square = optimize_diag && (csp.n_rows() == csp.n_cols());
+  std::vector<unsigned int> row_lengths (csp.n_rows());
+  for (unsigned int i=0; i<csp.n_rows(); ++i)
+    row_lengths[i] = csp.row_length(i) +
+                     (is_square ? 1 : 0);
+  reinit (csp.n_rows(), csp.n_cols(), row_lengths, is_square);
+
+                                  // now enter all the elements into
+                                  // the matrix. note that if the
+                                  // matrix is quadratic, then we
+                                  // already have the diagonal
+                                  // element preallocated
+  for (unsigned int row = 0; row<csp.n_rows(); ++row)
+    {
+      unsigned int *cols = &colnums[rowstart[row]] + (is_square ? 1 : 0);
+      CompressedSetSparsityPattern::CSSPIterator col_num = csp.row_begin (row);
+
+      for (; col_num != csp.row_end (row); ++col_num)
+       {
+         const unsigned int col = *col_num;
+         //      Assert (col < csp.n_cols(), ExcInvalidIndex(col,csp.n_cols()));
+         
+         if ((col!=row) || !is_square)
+           *cols++ = col;
+       }
+    }
+
+                                  // finally compress
+                                  // everything. this also sorts the
+                                  // entries within each row
+  compress ();
+}
+
+
+
+
 template <typename number>
 void SparsityPattern::copy_from (const FullMatrix<number> &matrix,
                                 const bool optimize_diag)

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.