]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add a preliminary, untested version of the chunk sparsity pattern (need to commit...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 27 Jul 2008 11:53:46 +0000 (11:53 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 27 Jul 2008 11:53:46 +0000 (11:53 +0000)
git-svn-id: https://svn.dealii.org/trunk@16451 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/lac/include/lac/chunk_sparsity_pattern.h b/deal.II/lac/include/lac/chunk_sparsity_pattern.h
new file mode 100644 (file)
index 0000000..d273465
--- /dev/null
@@ -0,0 +1,855 @@
+//---------------------------------------------------------------------------
+//    $Id: chunk_sparsity_pattern.h 15432 2007-11-03 03:08:43Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2008 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__chunk_sparsity_pattern_h
+#define __deal2__chunk_sparsity_pattern_h
+
+
+#include <base/config.h>
+#include <base/exceptions.h>
+#include <base/subscriptor.h>
+#include <base/vector_slice.h>
+
+#include <lac/sparsity_pattern.h>
+
+#include <vector>
+#include <iostream>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/*! @addtogroup Sparsity
+ *@{
+ */
+
+
+/**
+ * Structure representing the sparsity pattern of a sparse matrix.
+ *
+ * This class is an example of the "static" type of @ref Sparsity.
+ *
+ * It uses the compressed row storage (CSR) format to store data.
+ *
+ * @author Wolfgang Bangerth, Guido Kanschat and others
+ */
+class ChunkSparsityPattern : public Subscriptor
+{
+  public:
+                                    /**
+                                     * 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.
+                                     */
+    ChunkSparsityPattern ();
+    
+                                    /**
+                                     * 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
+                                     * (ChunkSparsityPattern());</tt>,
+                                     * with <tt>v</tt> a vector of
+                                     * ChunkSparsityPattern 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.
+                                     */
+    ChunkSparsityPattern (const ChunkSparsityPattern &);
+
+                                    /**
+                                     * Initialize a rectangular
+                                     * matrix.
+                                     *
+                                     * @arg m number of rows
+                                     * @arg n number of columns
+                                     * @arg max_per_row maximum
+                                     * number of nonzero entries per row
+                                     *
+                                     * @arg optimize_diagonal store
+                                     * diagonal entries first in row;
+                                     * see optimize_diagonal(). This
+                                     * takes effect for quadratic
+                                     * matrices only.
+                                     */
+    ChunkSparsityPattern (const unsigned int m,
+                         const unsigned int n,
+                         const unsigned int max_chunks_per_row,
+                         const unsigned int chunk_size,
+                         const bool optimize_diagonal = true);
+
+                                    /**
+                                     * Initialize a rectangular
+                                     * matrix.
+                                     *
+                                     * @arg m number of rows
+                                     * @arg n number of columns
+                                     *
+                                     * @arg row_lengths possible
+                                     * number of nonzero entries for
+                                     * each row.  This vector must
+                                     * have one entry for each row.
+                                     *
+                                     * @arg optimize_diagonal store
+                                     * diagonal entries first in row;
+                                     * see optimize_diagonal(). This
+                                     * takes effect for quadratic
+                                     * matrices only.
+                                     */
+    ChunkSparsityPattern (const unsigned int               m,
+                         const unsigned int               n,
+                         const std::vector<unsigned int>& row_lengths,
+                         const unsigned int chunk_size,
+                         const bool optimize_diagonal = true);
+    
+                                    /**
+                                     * Initialize a quadratic matrix
+                                     * of dimension <tt>n</tt> with
+                                     * at most <tt>max_per_row</tt>
+                                     * nonzero entries per row.
+                                     *
+                                     * This constructor automatically
+                                     * enables optimized storage of
+                                     * diagonal elements. To avoid
+                                     * this, use the constructor
+                                     * taking row and column numbers
+                                     * separately.
+                                     */
+    ChunkSparsityPattern (const unsigned int n,
+                         const unsigned int max_per_row,
+                         const unsigned int chunk_size);
+
+                                    /**
+                                     * Initialize a quadratic matrix.
+                                     *
+                                     * @arg m number of rows and columns
+                                     *
+                                     * @arg row_lengths possible
+                                     * number of nonzero entries for
+                                     * each row.  This vector must
+                                     * have one entry for each row.
+                                     *
+                                     * @arg optimize_diagonal store
+                                     * diagonal entries first in row;
+                                     * see optimize_diagonal().
+                                     */
+    ChunkSparsityPattern (const unsigned int               m,
+                         const std::vector<unsigned int>& row_lengths,
+                         const unsigned int               chunk_size,
+                         const bool optimize_diagonal = true);
+    
+                                    /**
+                                     * Destructor.
+                                     */
+    ~ChunkSparsityPattern ();
+
+                                    /**
+                                     * 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.
+                                     */
+    ChunkSparsityPattern & operator = (const ChunkSparsityPattern &);
+    
+                                    /**
+                                     * Reallocate memory and set up data
+                                     * structures for a new matrix with
+                                     * <tt>m </tt>rows and <tt>n</tt> columns,
+                                     * with at most <tt>max_per_row</tt>
+                                     * nonzero entries per row.
+                                     *
+                                     * This function simply maps its
+                                     * operations to the other
+                                     * <tt>reinit</tt> function.
+                                     */
+    void reinit (const unsigned int m,
+                const unsigned int n,
+                const unsigned int max_per_row,
+                const unsigned int chunk_size,
+                const bool optimize_diagonal = true);
+
+                                    /**
+                                     * Reallocate memory for a matrix
+                                     * of size <tt>m x n</tt>. The
+                                     * number of entries for each row
+                                     * is taken from the array
+                                     * <tt>row_lengths</tt> which has to
+                                     * give this number of each row
+                                     * <tt>i=1...m</tt>.
+                                     *
+                                     * If <tt>m*n==0</tt> all memory is freed,
+                                     * resulting in a total reinitialization
+                                     * of the object. If it is nonzero, new
+                                     * memory is only allocated if the new
+                                     * size extends the old one. This is done
+                                     * to save time and to avoid fragmentation
+                                     * of the heap.
+                                     *
+                                     * If the number of rows equals
+                                     * the number of columns and the
+                                     * last parameter is true,
+                                     * diagonal elements are stored
+                                     * first in each row to allow
+                                     * optimized access in relaxation
+                                     * methods of SparseMatrix.
+                                     */
+    void reinit (const unsigned int               m,
+                const unsigned int               n,
+                const std::vector<unsigned int> &row_lengths,
+                const unsigned int chunk_size,
+                const bool optimize_diagonal = true);
+
+                                    /**
+                                     * Same as above, but with a
+                                     * VectorSlice argument instead.
+                                     */
+    void reinit (const unsigned int               m,
+                const unsigned int               n,
+                const VectorSlice<const std::vector<unsigned int> > &row_lengths,
+                const unsigned int chunk_size,
+                const bool optimize_diagonal = true);
+    
+                                    /**
+                                     * This function compresses the sparsity
+                                     * structure that this object represents.
+                                     * It does so by eliminating unused
+                                     * entries and sorting the remaining ones
+                                     * to allow faster access by usage of
+                                     * binary search algorithms. A special
+                                     * sorting scheme is used for the
+                                     * diagonal entry of quadratic matrices,
+                                     * which is always the first entry of
+                                     * each row.
+                                     *
+                                     * The memory which is no more
+                                     * needed is released.
+                                     *
+                                     * SparseMatrix objects require the
+                                     * ChunkSparsityPattern objects they are
+                                     * initialized with to be compressed, to
+                                     * reduce memory requirements.
+                                     */
+    void compress ();
+    
+                                    /**
+                                     * This function can be used as a
+                                     * replacement for reinit(),
+                                     * subsequent calls to add() and
+                                     * a final call to close() if you
+                                     * know exactly in advance the
+                                     * entries that will form the
+                                     * matrix sparsity pattern.
+                                     *
+                                     * The first two parameters
+                                     * determine the size of the
+                                     * matrix. For the two last ones,
+                                     * note that a sparse matrix can
+                                     * be described by a sequence of
+                                     * rows, each of which is
+                                     * represented by a sequence of
+                                     * pairs of column indices and
+                                     * values. In the present
+                                     * context, the begin() and
+                                     * end() parameters designate
+                                     * iterators (of forward iterator
+                                     * type) into a container, one
+                                     * representing one row. The
+                                     * distance between begin()
+                                     * and end() should therefore
+                                     * be equal to
+                                     * n_rows(). These iterators
+                                     * may be iterators of
+                                     * <tt>std::vector</tt>,
+                                     * <tt>std::list</tt>, pointers into a
+                                     * C-style array, or any other
+                                     * iterator satisfying the
+                                     * requirements of a forward
+                                     * iterator. The objects pointed
+                                     * to by these iterators
+                                     * (i.e. what we get after
+                                     * applying <tt>operator*</tt> or
+                                     * <tt>operator-></tt> to one of these
+                                     * iterators) must be a container
+                                     * itself that provides functions
+                                     * <tt>begin</tt> and <tt>end</tt>
+                                     * designating a range of
+                                     * iterators that describe the
+                                     * contents of one
+                                     * line. Dereferencing these
+                                     * inner iterators must either
+                                     * yield a pair of an unsigned
+                                     * integer as column index and a
+                                     * value of arbitrary type (such
+                                     * a type would be used if we
+                                     * wanted to describe a sparse
+                                     * matrix with one such object),
+                                     * or simply an unsigned integer
+                                     * (of we only wanted to describe
+                                     * a sparsity pattern). The
+                                     * function is able to determine
+                                     * itself whether an unsigned
+                                     * integer or a pair is what we
+                                     * get after dereferencing the
+                                     * inner iterators, through some
+                                     * template magic.
+                                     *
+                                     * While the order of the outer
+                                     * iterators denotes the
+                                     * different rows of the matrix,
+                                     * the order of the inner
+                                     * iterator denoting the columns
+                                     * does not matter, as they are
+                                     * sorted internal to this
+                                     * function anyway.
+                                     *
+                                     * Since that all sounds very
+                                     * complicated, consider the
+                                     * following example code, which
+                                     * may be used to fill a sparsity
+                                     * pattern:
+                                     * @code
+                                     * std::vector<std::vector<unsigned int> > column_indices (n_rows);
+                                     * for (unsigned int row=0; row<n_rows; ++row)
+                                     *         // generate necessary columns in this row
+                                     *   fill_row (column_indices[row]);
+                                     *
+                                     * sparsity.copy_from (n_rows, n_cols,
+                                     *                     column_indices.begin(),
+                                     *                     column_indices.end());
+                                     * @endcode
+                                     *
+                                     * Note that this example works
+                                     * since the iterators
+                                     * dereferenced yield containers
+                                     * with functions <tt>begin</tt> and
+                                     * <tt>end</tt> (namely
+                                     * <tt>std::vector</tt>s), and the
+                                     * inner iterators dereferenced
+                                     * yield unsigned integers as
+                                     * column indices. Note that we
+                                     * could have replaced each of
+                                     * the two <tt>std::vector</tt>
+                                     * occurrences by <tt>std::list</tt>,
+                                     * and the inner one by
+                                     * <tt>std::set</tt> as well.
+                                     *
+                                     * Another example would be as
+                                     * follows, where we initialize a
+                                     * whole matrix, not only a
+                                     * sparsity pattern:
+                                     * @code
+                                     * std::vector<std::map<unsigned int,double> > entries (n_rows);
+                                     * for (unsigned int row=0; row<n_rows; ++row)
+                                     *         // generate necessary pairs of columns
+                                     *         // and corresponding values in this row
+                                     *   fill_row (entries[row]);
+                                     *
+                                     * sparsity.copy_from (n_rows, n_cols,
+                                     *                     column_indices.begin(),
+                                     *                     column_indices.end());
+                                     * matrix.reinit (sparsity);
+                                     * matrix.copy_from (column_indices.begin(),
+                                     *                   column_indices.end());
+                                     * @endcode
+                                     *
+                                     * This example works because
+                                     * dereferencing iterators of the
+                                     * inner type yields a pair of
+                                     * unsigned integers and a value,
+                                     * the first of which we take as
+                                     * column index. As previously,
+                                     * the outer <tt>std::vector</tt>
+                                     * could be replaced by
+                                     * <tt>std::list</tt>, and the inner
+                                     * <tt>std::map<unsigned int,double></tt>
+                                     * could be replaced by
+                                     * <tt>std::vector<std::pair<unsigned int,double> ></tt>,
+                                     * or a list or set of such
+                                     * pairs, as they all return
+                                     * iterators that point to such
+                                     * pairs.
+                                     */
+    template <typename ForwardIterator>
+    void copy_from (const unsigned int    n_rows,
+                   const unsigned int    n_cols,
+                   const ForwardIterator begin,
+                   const ForwardIterator end,
+                   const unsigned int chunk_size,
+                   const bool optimize_diagonal = true);
+
+                                    /**
+                                     * Copy data from an object of
+                                     * type
+                                     * CompressedSparsityPattern.
+                                     * Previous content of this
+                                     * object is lost, and the
+                                     * sparsity pattern is in
+                                     * compressed mode afterwards.
+                                     */
+    void copy_from (const CompressedSparsityPattern &csp,
+                   const unsigned int chunk_size,
+                   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 unsigned int chunk_size,
+                   const bool optimize_diagonal = true);
+
+
+                                    /**
+                                     * Take a full matrix and use its
+                                     * nonzero entries to generate a
+                                     * sparse matrix entry pattern
+                                     * for this object.
+                                     *
+                                     * Previous content of this
+                                     * object is lost, and the
+                                     * sparsity pattern is in
+                                     * compressed mode afterwards.
+                                     */
+    template <typename number>
+    void copy_from (const FullMatrix<number> &matrix,
+                   const unsigned int chunk_size,
+                   const bool optimize_diagonal = true);
+    
+                                    /**
+                                     * 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. Before compression, this equals the
+                                     * number given to the constructor, while
+                                     * after compression, it equals the maximum
+                                     * number of entries actually allocated by
+                                     * the user.
+                                     */
+    unsigned int max_entries_per_row () const;
+
+                                    /**
+                                     * Add a nonzero entry to the matrix.
+                                     * This function may only be called
+                                     * for non-compressed sparsity patterns.
+                                     *
+                                     * If the entry already exists, nothing
+                                     * bad happens.
+                                     */
+    void add (const unsigned int i, 
+             const unsigned int j);
+    
+                                    /**
+                                     * 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
+                                     * quadratic matrix.
+                                     */
+    void symmetrize ();
+    
+                                    /**
+                                     * Return number of rows of this
+                                     * matrix, which equals the dimension
+                                     * of the image space.
+                                     */
+    inline unsigned int n_rows () const;
+
+                                    /**
+                                     * Return number of columns of this
+                                     * matrix, which equals the dimension
+                                     * of the range space.
+                                     */
+    inline unsigned int n_cols () const;
+
+                                    /**
+                                     * Check if a value at a certain
+                                     * position may be non-zero.
+                                     */
+    bool exists (const unsigned int i,
+                 const unsigned int j) 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. Consequently, the maximum
+                                     * bandwidth a $n\times m$ matrix can
+                                     * have is $\max\{n-1,m-1\}$.
+                                     */
+    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;
+
+                                    /**
+                                     * Return whether the structure is 
+                                     * compressed or not.
+                                     */
+    bool is_compressed () const;
+
+                                    /**
+                                     * Determine whether the matrix
+                                     * uses special convention for
+                                     * quadratic matrices.
+                                     *
+                                     * A return value <tt>true</tt> means
+                                     * that diagonal elements are stored
+                                     * first in each row. A number of
+                                     * functions in this class and the
+                                     * library in general, for example
+                                     * relaxation methods like Jacobi() and
+                                     * SOR(), require this to make their
+                                     * operations more efficient, since they
+                                     * need to quickly access the diagonal
+                                     * elements and do not have to search for
+                                     * them if they are the first element of
+                                     * each row. A side effect of this scheme
+                                     * is that each row contains at least one
+                                     * element, even if the row is empty
+                                     * (i.e. the diagonal element exists, but
+                                     * has value zero).
+                                     *
+                                     * A return value <tt>false</tt> means
+                                     * that diagonal elements are stored
+                                     * anywhere in the row, or not at all. In
+                                     * particular, a row or even the whole
+                                     * matrix may be empty. This can be used
+                                     * if you have block matrices where the
+                                     * off-diagonal blocks are quadratic but
+                                     * are never used for operations like the
+                                     * ones mentioned above. In this case,
+                                     * some memory can be saved by not using
+                                     * the diagonal storage optimization.
+                                      */
+    bool optimize_diagonal () const;
+
+    
+                                    /**
+                                     * Write the data of this object
+                                     * en bloc to a file. This is
+                                     * done in a binary mode, so the
+                                     * output is neither readable by
+                                     * humans nor (probably) by other
+                                     * computers using a different
+                                     * operating system of number
+                                     * format.
+                                     *
+                                     * The purpose of this function
+                                     * is that you can swap out
+                                     * matrices and sparsity pattern
+                                     * if you are short of memory,
+                                     * want to communicate between
+                                     * different programs, or allow
+                                     * objects to be persistent
+                                     * across different runs of the
+                                     * program.
+                                     */
+    void block_write (std::ostream &out) const;
+
+                                    /**
+                                     * Read data that has previously
+                                     * been written by block_write()
+                                     * from a file. This is done
+                                     * using the inverse operations
+                                     * to the above function, so it
+                                     * is reasonably fast because the
+                                     * bitstream is not interpreted
+                                     * except for a few numbers up
+                                     * front.
+                                     *
+                                     * The object is resized on this
+                                     * operation, and all previous
+                                     * contents are lost.
+                                     *
+                                     * A primitive form of error
+                                     * checking is performed which
+                                     * will recognize the bluntest
+                                     * attempts to interpret some
+                                     * data as a vector stored
+                                     * bitwise to a file, but not
+                                     * more.
+                                     */
+    void block_read (std::istream &in);
+
+                                    /**
+                                     * Print the sparsity of the
+                                     * matrix. The output consists of
+                                     * one line per row of the format
+                                     * <tt>[i,j1,j2,j3,...]</tt>. <i>i</i>
+                                     * is the row number and
+                                     * <i>jn</i> are the allocated
+                                     * columns in this row.
+                                     */
+    void print (std::ostream &out) const;
+
+                                    /**
+                                     * Print the sparsity of the matrix
+                                     * in a format that <tt>gnuplot</tt> 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
+                                     * <tt>plot</tt> command.
+                                     */
+    void print_gnuplot (std::ostream &out) const;
+
+                                    /**
+                                     * Determine an estimate for the
+                                     * memory consumption (in bytes)
+                                     * of this object. See
+                                     * MemoryConsumption.
+                                     */
+    unsigned int memory_consumption () const;
+
+                                    /** @addtogroup Exceptions
+                                     * @{ */
+                                    /**
+                                     * Exception
+                                     */
+    DeclException1 (ExcInvalidNumber,
+                   int,
+                   << "The provided number is invalid here: " << arg1);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidIndex,
+                   int, int,
+                   << "The given index " << arg1
+                   << " should be less than " << arg2 << ".");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcNotEnoughSpace,
+                   int, int,
+                   << "Upon entering a new entry to row " << arg1
+                   << ": there was no free entry any more. " << std::endl
+                   << "(Maximum number of entries for this row: "
+                   << arg2 << "; maybe the matrix is already compressed?)");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNotCompressed);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcMatrixIsCompressed);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcEmptyObject);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInvalidConstructorCall);
+                                    /**
+                                     * This exception is thrown if
+                                     * the matrix does not follow the
+                                     * convention of storing diagonal
+                                     * elements first in row. Refer
+                                     * to
+                                     * SparityPattern::optimize_diagonal()
+                                     * for more information.
+                                     */
+    DeclException0 (ExcDiagonalNotOptimized);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcIteratorRange,
+                   int, int,
+                   << "The iterators denote a range of " << arg1
+                   << " elements, but the given number of rows was " << arg2);
+                                     /**
+                                      * Exception
+                                      */
+    DeclException0 (ExcMETISNotInstalled);
+                                     /**
+                                      * Exception
+                                      */
+    DeclException1 (ExcInvalidNumberOfPartitions,
+                    int,
+                    << "The number of partitions you gave is " << arg1
+                    << ", but must be greater than zero.");
+                                     /**
+                                      * Exception
+                                      */
+    DeclException2 (ExcInvalidArraySize,
+                    int, int,
+                    << "The array has size " << arg1 << " but should have size "
+                    << arg2);
+                                    //@}
+  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;
+
+                                    /**
+                                     * The size of chunks.
+                                     */
+    unsigned int chunk_size;
+
+                                    /**
+                                     * The reduced sparsity pattern. We store
+                                     * only which chunks exist, with each
+                                     * chunk a block in the matrix of size
+                                     * chunk_size by chunk_size.
+                                     */
+    SparsityPattern sparsity_pattern;
+};
+
+
+/*@}*/
+/*---------------------- Inline functions -----------------------------------*/
+
+#ifndef DOXYGEN
+
+
+inline
+unsigned int
+ChunkSparsityPattern::n_rows () const
+{
+  return rows;
+}
+
+
+inline
+unsigned int
+ChunkSparsityPattern::n_cols () const
+{
+  return cols;
+}
+
+
+inline
+bool
+ChunkSparsityPattern::is_compressed () const
+{
+  return sparsity_pattern.compressed;
+}
+
+
+inline
+bool
+ChunkSparsityPattern::optimize_diagonal () const
+{
+  return sparsity_pattern.diagonal_optimized;
+}
+
+
+inline
+unsigned int
+ChunkSparsityPattern::n_nonzero_elements () const
+{
+  return (sparsity_pattern.n_nonzero_elements() *
+         chunk_size *
+         chunk_size);
+}
+
+
+
+template <typename ForwardIterator>
+void
+ChunkSparsityPattern::copy_from (const unsigned int    n_rows,
+                                const unsigned int    n_cols,
+                                const ForwardIterator begin,
+                                const ForwardIterator end,
+                                const unsigned int    chunk_size,
+                                const bool optimize_diag)
+{
+  Assert (static_cast<unsigned int>(std::distance (begin, end)) == n_rows,
+         ExcIteratorRange (std::distance (begin, end), n_rows));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+#endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 1378a3a3d0c2a8ae2011aad3f55788554fa29f8c..2209cdf516f703604d8dee2f3c9b24616d506979 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -24,6 +24,7 @@
 DEAL_II_NAMESPACE_OPEN
 
 class SparsityPattern;
+class ChunkSparsityPattern;
 template <typename number> class FullMatrix;
 template <typename number> class SparseMatrix;
 template <class VECTOR> class VectorSlice;
@@ -148,6 +149,42 @@ namespace internals
            len = half;
        }
     }
+
+
+                                    /**
+                                     * Helper function to get the
+                                     * column index from a
+                                     * dereferenced iterator in the
+                                     * copy_from() function, if
+                                     * the inner iterator type points
+                                     * to plain unsigned integers.
+                                     */
+    unsigned int
+    get_column_index_from_iterator (const unsigned int i);
+    
+                                    /**
+                                     * Helper function to get the
+                                     * column index from a
+                                     * dereferenced iterator in the
+                                     * copy_from() function, if
+                                     * the inner iterator type points
+                                     * to pairs of unsigned integers
+                                     * and some other value.
+                                     */
+    template <typename value>
+    unsigned int
+    get_column_index_from_iterator (const std::pair<unsigned int, value> &i);
+
+                                    /**
+                                     * Likewise, but sometimes needed
+                                     * for certain types of
+                                     * containers that make the first
+                                     * element of the pair constant
+                                     * (such as <tt>std::map</tt>).
+                                     */
+    template <typename value>
+    unsigned int
+    get_column_index_from_iterator (const std::pair<const unsigned int, value> &i);
     
   }
   
@@ -1069,10 +1106,12 @@ class SparsityPattern : public Subscriptor
                                     /**
                                      * 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.
+                                     * bandwidth is the maximum of $|i-j|$
+                                     * for which the index pair $(i,j)$
+                                     * represents a nonzero entry of the
+                                     * matrix. Consequently, the maximum
+                                     * bandwidth a $n\times m$ matrix can
+                                     * have is $\max\{n-1,m-1\}$.
                                      */
     unsigned int bandwidth () const;
 
@@ -1556,49 +1595,15 @@ class SparsityPattern : public Subscriptor
                                      * Is special treatment of
                                      * diagonals enabled?
                                      */
-    bool diagonal_optimized;
-    
-                                    /**
-                                     * Helper function to get the
-                                     * column index from a
-                                     * dereferenced iterator in the
-                                     * copy_from() function, if
-                                     * the inner iterator type points
-                                     * to plain unsigned integers.
-                                     */
-    static
-    unsigned int
-    get_column_index_from_iterator (const unsigned int i);
-    
-                                    /**
-                                     * Helper function to get the
-                                     * column index from a
-                                     * dereferenced iterator in the
-                                     * copy_from() function, if
-                                     * the inner iterator type points
-                                     * to pairs of unsigned integers
-                                     * and some other value.
-                                     */
-    template <typename value>
-    unsigned int
-    get_column_index_from_iterator (const std::pair<unsigned int, value> &i);
-
-                                    /**
-                                     * Likewise, but sometimes needed
-                                     * for certain types of
-                                     * containers that make the first
-                                     * element of the pair constant
-                                     * (such as <tt>std::map</tt>).
-                                     */
-    template <typename value>
-    unsigned int
-    get_column_index_from_iterator (const std::pair<const unsigned int, value> &i);
+    bool diagonal_optimized;    
     
                                     /**
                                      * Make all sparse matrices
                                      * friends of this class.
                                      */
     template <typename number> friend class SparseMatrix;
+
+    friend class ChunkSparsityPattern;
 };
 
 
@@ -1946,34 +1951,37 @@ SparsityPattern::n_nonzero_elements () const
 
 
 
-inline
-unsigned int
-SparsityPattern::
-get_column_index_from_iterator (const unsigned int i)
+namespace internal
 {
-  return i;
-}
+  namespace SparsityPatternTools
+  {
+    inline
+    unsigned int
+    get_column_index_from_iterator (const unsigned int i)
+    {
+      return i;
+    }
 
 
 
-template <typename value>
-inline
-unsigned int
-SparsityPattern::
-get_column_index_from_iterator (const std::pair<unsigned int, value> &i)
-{
-  return i.first;
-}
+    template <typename value>
+    inline
+    unsigned int
+    get_column_index_from_iterator (const std::pair<unsigned int, value> &i)
+    {
+      return i.first;
+    }
 
 
 
-template <typename value>
-inline
-unsigned int
-SparsityPattern::
-get_column_index_from_iterator (const std::pair<const unsigned int, value> &i)
-{
-  return i.first;
+    template <typename value>
+    inline
+    unsigned int
+    get_column_index_from_iterator (const std::pair<const unsigned int, value> &i)
+    {
+      return i.first;
+    }
+  }
 }
 
 
@@ -2025,7 +2033,8 @@ SparsityPattern::copy_from (const unsigned int    n_rows,
       const inner_iterator end_of_row = i->end();
       for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
        {
-         const unsigned int col = get_column_index_from_iterator(*j);
+         const unsigned int col
+           = internal::SparsityPatternTools::get_column_index_from_iterator(*j);
          Assert (col < n_cols, ExcInvalidIndex(col,n_cols));
          
          if ((col!=row) || !is_square)
diff --git a/deal.II/lac/source/chunk_sparsity_pattern.cc b/deal.II/lac/source/chunk_sparsity_pattern.cc
new file mode 100644 (file)
index 0000000..b1e3216
--- /dev/null
@@ -0,0 +1,438 @@
+//---------------------------------------------------------------------------
+//    $Id: chunk_sparsity_pattern.cc 14783 2007-06-18 14:52:01Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2008 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/chunk_sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/compressed_set_sparsity_pattern.h>
+#include <lac/full_matrix.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+ChunkSparsityPattern::ChunkSparsityPattern ()
+{
+  reinit (0,0,0,0);
+}
+
+
+
+ChunkSparsityPattern::ChunkSparsityPattern (const ChunkSparsityPattern &s)
+                :
+               Subscriptor(),
+               chunk_size (s.chunk_size),
+               sparsity_pattern(s.sparsity_pattern)
+{
+  Assert (s.rows == 0, ExcInvalidConstructorCall());
+  Assert (s.cols == 0, ExcInvalidConstructorCall());
+  
+  reinit (0,0,0,0, false);
+}
+
+
+
+ChunkSparsityPattern::ChunkSparsityPattern (const unsigned int m,
+                                           const unsigned int n,
+                                           const unsigned int max_per_row,
+                                           const unsigned int chunk_size,
+                                           const bool optimize_diag)
+{
+  reinit (m,n,max_per_row, chunk_size, optimize_diag);
+}
+
+
+
+ChunkSparsityPattern::ChunkSparsityPattern (
+  const unsigned int m,
+  const unsigned int n,
+  const std::vector<unsigned int>& row_lengths,
+  const unsigned int chunk_size,
+  const bool optimize_diag) 
+{
+  reinit (m, n, row_lengths, chunk_size, optimize_diag);
+}
+
+
+
+ChunkSparsityPattern::ChunkSparsityPattern (const unsigned int n,
+                                           const unsigned int max_per_row,
+                                           const unsigned int chunk_size)
+{
+  reinit (n, n, max_per_row, chunk_size, true);
+}
+
+
+
+ChunkSparsityPattern::ChunkSparsityPattern (
+  const unsigned int               m,
+  const std::vector<unsigned int>& row_lengths,
+  const unsigned int chunk_size,
+  const bool optimize_diag) 
+{
+  reinit (m, m, row_lengths, chunk_size, optimize_diag);
+}
+
+
+
+ChunkSparsityPattern::~ChunkSparsityPattern ()
+{}
+
+
+
+ChunkSparsityPattern &
+ChunkSparsityPattern::operator = (const ChunkSparsityPattern &s)
+{
+  Assert (s.rows == 0, ExcInvalidConstructorCall());
+  Assert (s.cols == 0, ExcInvalidConstructorCall());
+
+                                  // perform the checks in the underlying
+                                  // object as well
+  sparsity_pattern = s.sparsity_pattern;
+  
+  return *this;
+}
+
+
+
+void
+ChunkSparsityPattern::reinit (const unsigned int m,
+                             const unsigned int n,
+                             const unsigned int max_per_row,
+                             const unsigned int chunk_size,
+                             const bool optimize_diag)
+{
+                                  // simply map this function to the
+                                  // other @p{reinit} function
+  const std::vector<unsigned int> row_lengths (m, max_per_row);
+  reinit (m, n, row_lengths, chunk_size, optimize_diag);
+}
+
+
+
+void
+ChunkSparsityPattern::reinit (
+  const unsigned int m,
+  const unsigned int n,
+  const VectorSlice<const std::vector<unsigned int> >&row_lengths,
+  const unsigned int chunk_size,
+  const bool optimize_diag)
+{
+  Assert (row_lengths.size() == m, ExcInvalidNumber (m));
+         
+  rows = m;
+  cols = n;
+
+  this->chunk_size = chunk_size;
+
+                                  // pass down to the necessary information
+                                  // to the underlying object. we need to
+                                  // calculate how many chunks we need: we
+                                  // need to round up (m/chunk_size) and
+                                  // (n/chunk_size). rounding up in integer
+                                  // arithmetic equals
+                                  // ((m+chunk_size-1)/chunk_size):
+  const unsigned int m_chunks = (m+chunk_size) / chunk_size,
+                    n_chunks = (n+chunk_size) / chunk_size;  
+
+                                  // compute the maximum number of chunks in
+                                  // each row. the passed array denotes the
+                                  // number of entries in each row -- in the
+                                  // worst case, these are all in independent
+                                  // chunks, so we have to calculate it as
+                                  // follows:
+  std::vector<unsigned int> chunk_row_lengths (m_chunks, 0);
+  for (unsigned int i=0; i<m; ++i)
+    chunk_row_lengths[i/chunk_size]
+      = std::max (chunk_row_lengths[i/chunk_size],
+                 row_lengths[i]);
+
+  sparsity_pattern.reinit (m_chunks,
+                          n_chunks,
+                          chunk_row_lengths,
+                          optimize_diag);
+}
+
+
+
+void
+ChunkSparsityPattern::compress ()
+{
+  sparsity_pattern.compress ();
+}
+
+
+
+void
+ChunkSparsityPattern::copy_from (const CompressedSparsityPattern &csp,
+                                const unsigned int chunk_size,
+                                const bool         optimize_diag) 
+{
+                                  // count number of entries per row, then
+                                  // initialize the underlying sparsity
+                                  // pattern
+  std::vector<unsigned int> entries_per_row (csp.n_rows(), 0);
+  for (unsigned int row = 0; row<csp.n_rows(); ++row)
+    for (unsigned int j=0; j<csp.row_length(row); ++j)
+      ++entries_per_row[row];
+
+  reinit (csp.n_rows(), csp.n_cols(),
+         entries_per_row,
+         chunk_size, optimize_diag);
+      
+                                  // then actually fill it  
+  for (unsigned int row = 0; row<csp.n_rows(); ++row)
+    for (unsigned int j=0; j<csp.row_length(row); ++j)
+      add (row, csp.column_number(row,j));
+
+                                  // finally compress
+  compress ();
+}
+
+
+
+void
+ChunkSparsityPattern::copy_from (const CompressedSetSparsityPattern &csp,
+                                const unsigned int chunk_size,
+                                const bool         optimize_diag) 
+{
+                                  // count number of entries per row, then
+                                  // initialize the underlying sparsity
+                                  // pattern
+  std::vector<unsigned int> entries_per_row (csp.n_rows(), 0);
+  for (unsigned int row = 0; row<csp.n_rows(); ++row)
+    {
+      CompressedSetSparsityPattern::row_iterator col_num = csp.row_begin (row);
+
+      for (; col_num != csp.row_end (row); ++col_num)
+       ++entries_per_row[row];
+    }
+
+  reinit (csp.n_rows(), csp.n_cols(),
+         entries_per_row,
+         chunk_size, optimize_diag);
+      
+                                  // then actually fill it  
+  for (unsigned int row = 0; row<csp.n_rows(); ++row)
+    {
+      CompressedSetSparsityPattern::row_iterator col_num = csp.row_begin (row);
+
+      for (; col_num != csp.row_end (row); ++col_num)
+       add (row, *col_num);
+    }
+  
+                                  // finally compress
+  compress ();
+}
+
+
+
+
+template <typename number>
+void ChunkSparsityPattern::copy_from (const FullMatrix<number> &matrix,
+                                     const unsigned int chunk_size,
+                                     const bool         optimize_diag)
+{
+                                  // count number of entries per row, then
+                                  // initialize the underlying sparsity
+                                  // pattern
+  std::vector<unsigned int> entries_per_row (matrix.m(), 0);
+  for (unsigned int row=0; row<matrix.m(); ++row)
+    for (unsigned int col=0; col<matrix.n(); ++col)
+      if (matrix(row,col) != 0)
+       ++entries_per_row[row];
+
+  reinit (matrix.m(), matrix.n(),
+         entries_per_row,
+         chunk_size, optimize_diag);
+      
+                                  // then actually fill it  
+  for (unsigned int row=0; row<matrix.m(); ++row)
+    for (unsigned int col=0; col<matrix.n(); ++col)
+      if (matrix(row,col) != 0)
+       add (row,col);
+
+                                  // finally compress
+  compress ();
+}
+
+
+void
+ChunkSparsityPattern::reinit (
+  const unsigned int m,
+  const unsigned int n,
+  const std::vector<unsigned int>& row_lengths,
+  const unsigned int chunk_size,
+  const bool optimize_diag)
+{
+  reinit(m, n, make_slice(row_lengths), chunk_size, optimize_diag);
+}
+
+
+
+bool
+ChunkSparsityPattern::empty () const
+{
+  return sparsity_pattern.empty();
+}
+
+
+
+unsigned int
+ChunkSparsityPattern::max_entries_per_row () const 
+{
+  return sparsity_pattern.max_entries_per_row() * chunk_size;
+}
+
+
+
+void
+ChunkSparsityPattern::add (const unsigned int i,
+                          const unsigned int j)
+{
+  Assert (i<rows, ExcInvalidIndex(i,rows));
+  Assert (j<cols, ExcInvalidIndex(j,cols));
+
+  sparsity_pattern.add (i/chunk_size, j/chunk_size);
+}
+
+
+bool
+ChunkSparsityPattern::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 sparsity_pattern.exists (i/chunk_size,
+                                 j/chunk_size);
+}
+
+
+
+void
+ChunkSparsityPattern::symmetrize () 
+{
+                                  // matrix must be square. note that the for
+                                  // some matrix sizes, the current sparsity
+                                  // pattern may not be square even if the
+                                  // underlying sparsity pattern is (e.g. a
+                                  // 10x11 matrix with chunk_size 4)
+  Assert (rows==cols, ExcNotQuadratic());
+
+  sparsity_pattern.symmetrize ();
+}
+
+
+
+void
+ChunkSparsityPattern::print (std::ostream &/*out*/) const
+{
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+void
+ChunkSparsityPattern::print_gnuplot (std::ostream &/*out*/) const
+{
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+unsigned int
+ChunkSparsityPattern::bandwidth () const
+{
+                                  // calculate the bandwidth from that of the
+                                  // underlying sparsity pattern. note that
+                                  // even if the bandwidth of that is zero,
+                                  // then the bandwidth of the chunky pattern
+                                  // is chunk_size-1, if it is 1 then the
+                                  // chunky pattern has
+                                  // chunk_size+(chunk_size-1), etc
+                                  //
+                                  // we'll cut it off at max(n(),m())
+  return std::min (sparsity_pattern.bandwidth()*chunk_size
+                  + (chunk_size-1),
+                  std::max(n_rows(), n_cols()));
+}
+
+
+void
+ChunkSparsityPattern::block_write (std::ostream &out) const 
+{
+  AssertThrow (out, ExcIO());
+
+                                   // first the simple objects,
+                                   // bracketed in [...]
+  out << '['
+      << rows << ' '
+      << cols << ' '
+      << "][";
+                                  // then the underlying sparsity pattern
+  sparsity_pattern.block_write (out);
+  out << ']';
+  
+  AssertThrow (out, ExcIO());  
+}
+
+
+
+void
+ChunkSparsityPattern::block_read (std::istream &in)
+{
+  AssertThrow (in, ExcIO());
+
+  char c;
+
+                                   // first read in simple data
+  in >> c;
+  AssertThrow (c == '[', ExcIO());
+  in >> rows
+     >> cols;
+
+  in >> c;
+  AssertThrow (c == ']', ExcIO());
+  in >> c;
+  AssertThrow (c == '[', ExcIO());
+
+                                  // then read the underlying sparsity
+                                  // pattern
+  sparsity_pattern.block_read (in);
+
+  in >> c;
+  AssertThrow (c == ']', ExcIO());  
+}
+
+
+
+unsigned int
+ChunkSparsityPattern::memory_consumption () const
+{
+  return (sizeof(*this) +
+         sparsity_pattern.memory_consumption());
+}
+
+
+
+// explicit instantiations
+template
+void ChunkSparsityPattern::copy_from<float> (const FullMatrix<float> &,
+                                            const unsigned int,
+                                            bool);
+template
+void ChunkSparsityPattern::copy_from<double> (const FullMatrix<double> &,
+                                             const unsigned int,
+                                             bool);
+
+DEAL_II_NAMESPACE_CLOSE

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.