From: bangerth Date: Sun, 27 Jul 2008 11:53:46 +0000 (+0000) Subject: Add a preliminary, untested version of the chunk sparsity pattern (need to commit... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d6feb84b698748d8bf2fa6e63a43d06ab0342140;p=dealii-svn.git Add a preliminary, untested version of the chunk sparsity pattern (need to commit because the laptop this is on will have its harddrive erased in a couple of hours). git-svn-id: https://svn.dealii.org/trunk@16451 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 index 0000000000..d2734656d2 --- /dev/null +++ b/deal.II/lac/include/lac/chunk_sparsity_pattern.h @@ -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 +#include +#include +#include + +#include + +#include +#include + +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 + * v.push_back + * (ChunkSparsityPattern());, + * with v 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& row_lengths, + const unsigned int chunk_size, + const bool optimize_diagonal = true); + + /** + * Initialize a quadratic matrix + * of dimension n with + * at most max_per_row + * 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& 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 + * m rows and n columns, + * with at most max_per_row + * nonzero entries per row. + * + * This function simply maps its + * operations to the other + * reinit 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 m x n. The + * number of entries for each row + * is taken from the array + * row_lengths which has to + * give this number of each row + * i=1...m. + * + * If m*n==0 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 &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 > &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 + * std::vector, + * std::list, 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 operator* or + * operator-> to one of these + * iterators) must be a container + * itself that provides functions + * begin and end + * 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 > column_indices (n_rows); + * for (unsigned int row=0; rowbegin and + * end (namely + * std::vectors), and the + * inner iterators dereferenced + * yield unsigned integers as + * column indices. Note that we + * could have replaced each of + * the two std::vector + * occurrences by std::list, + * and the inner one by + * std::set as well. + * + * Another example would be as + * follows, where we initialize a + * whole matrix, not only a + * sparsity pattern: + * @code + * std::vector > entries (n_rows); + * for (unsigned int row=0; rowstd::vector + * could be replaced by + * std::list, and the inner + * std::map + * could be replaced by + * std::vector >, + * or a list or set of such + * pairs, as they all return + * iterators that point to such + * pairs. + */ + template + 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 + void copy_from (const FullMatrix &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 true 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 false 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 + * [i,j1,j2,j3,...]. i + * is the row number and + * jn are the allocated + * columns in this row. + */ + void print (std::ostream &out) const; + + /** + * Print the sparsity of the matrix + * in a format that gnuplot understands + * and which can be used to plot the + * sparsity pattern in a graphical + * way. The format consists of pairs + * i j 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 + * (0,0) 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 + * plot 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 +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(std::distance (begin, end)) == n_rows, + ExcIteratorRange (std::distance (begin, end), n_rows)); + + Assert (false, ExcNotImplemented()); +} + + +#endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/lac/include/lac/sparsity_pattern.h b/deal.II/lac/include/lac/sparsity_pattern.h index 1378a3a3d0..2209cdf516 100644 --- a/deal.II/lac/include/lac/sparsity_pattern.h +++ b/deal.II/lac/include/lac/sparsity_pattern.h @@ -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 class FullMatrix; template class SparseMatrix; template 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 + unsigned int + get_column_index_from_iterator (const std::pair &i); + + /** + * Likewise, but sometimes needed + * for certain types of + * containers that make the first + * element of the pair constant + * (such as std::map). + */ + template + unsigned int + get_column_index_from_iterator (const std::pair &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 - unsigned int - get_column_index_from_iterator (const std::pair &i); - - /** - * Likewise, but sometimes needed - * for certain types of - * containers that make the first - * element of the pair constant - * (such as std::map). - */ - template - unsigned int - get_column_index_from_iterator (const std::pair &i); + bool diagonal_optimized; /** * Make all sparse matrices * friends of this class. */ template 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 -inline -unsigned int -SparsityPattern:: -get_column_index_from_iterator (const std::pair &i) -{ - return i.first; -} + template + inline + unsigned int + get_column_index_from_iterator (const std::pair &i) + { + return i.first; + } -template -inline -unsigned int -SparsityPattern:: -get_column_index_from_iterator (const std::pair &i) -{ - return i.first; + template + inline + unsigned int + get_column_index_from_iterator (const std::pair &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 index 0000000000..b1e3216729 --- /dev/null +++ b/deal.II/lac/source/chunk_sparsity_pattern.cc @@ -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 +#include +#include +#include + + +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& 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& 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 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 >&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 chunk_row_lengths (m_chunks, 0); + for (unsigned int i=0; i entries_per_row (csp.n_rows(), 0); + for (unsigned int row = 0; row entries_per_row (csp.n_rows(), 0); + for (unsigned int row = 0; row +void ChunkSparsityPattern::copy_from (const FullMatrix &matrix, + const unsigned int chunk_size, + const bool optimize_diag) +{ + // count number of entries per row, then + // initialize the underlying sparsity + // pattern + std::vector entries_per_row (matrix.m(), 0); + for (unsigned int row=0; row& 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> 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 (const FullMatrix &, + const unsigned int, + bool); +template +void ChunkSparsityPattern::copy_from (const FullMatrix &, + const unsigned int, + bool); + +DEAL_II_NAMESPACE_CLOSE