--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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