include/deal.II/lac/block_sparsity_pattern.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: block_sparsity_pattern.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__block_sparsity_pattern_h
00013 #define __deal2__block_sparsity_pattern_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/table.h>
00019 #include <deal.II/base/subscriptor.h>
00020 #include <deal.II/base/smartpointer.h>
00021 #include <deal.II/lac/sparsity_pattern.h>
00022 #include <deal.II/lac/trilinos_sparsity_pattern.h>
00023 #include <deal.II/lac/compressed_sparsity_pattern.h>
00024 #include <deal.II/lac/compressed_set_sparsity_pattern.h>
00025 #include <deal.II/lac/compressed_simple_sparsity_pattern.h>
00026 #include <deal.II/lac/block_indices.h>
00027 
00028 DEAL_II_NAMESPACE_OPEN
00029 
00030 
00031 template <typename number> class BlockSparseMatrix;
00032 class BlockSparsityPattern;
00033 class BlockCompressedSparsityPattern;
00034 class BlockCompressedSimpleSparsityPattern;
00035 class BlockCompressedSetSparsityPattern;
00036 #ifdef DEAL_II_USE_TRILINOS
00037 namespace TrilinosWrappers
00038 {
00039   class BlockSparsityPattern;
00040 }
00041 #endif
00042 
00081 template <class SparsityPatternBase>
00082 class BlockSparsityPatternBase : public Subscriptor
00083 {
00084   public:
00097     static const unsigned int invalid_entry = SparsityPattern::invalid_entry;
00098 
00109     BlockSparsityPatternBase ();
00110 
00119     BlockSparsityPatternBase (const unsigned int n_block_rows,
00120                               const unsigned int n_block_columns);
00121 
00133     BlockSparsityPatternBase (const BlockSparsityPatternBase &bsp);
00134 
00138     ~BlockSparsityPatternBase ();
00139 
00166     void reinit (const unsigned int n_block_rows,
00167                  const unsigned int n_block_columns);
00168 
00177     BlockSparsityPatternBase & operator = (const BlockSparsityPatternBase &);
00178 
00191     void collect_sizes ();
00192 
00197     SparsityPatternBase &
00198     block (const unsigned int row,
00199            const unsigned int column);
00200 
00201 
00207     const SparsityPatternBase &
00208     block (const unsigned int row,
00209            const unsigned int column) const;
00210 
00217     const BlockIndices &
00218     get_row_indices () const;
00219 
00226     const BlockIndices &
00227     get_column_indices () const;
00228 
00236     void compress ();
00237 
00242     unsigned int n_block_rows () const;
00243 
00248     unsigned int n_block_cols () const;
00249 
00261     bool empty () const;
00262 
00271     unsigned int max_entries_per_row () const;
00272 
00285     void add (const unsigned int i, const unsigned int j);
00286 
00302     template <typename ForwardIterator>
00303     void add_entries (const unsigned int row,
00304                       ForwardIterator    begin,
00305                       ForwardIterator    end,
00306                       const bool         indices_are_sorted = false);
00307 
00316     unsigned int n_rows () const;
00317 
00326     unsigned int n_cols () const;
00327 
00332     bool exists (const unsigned int i, const unsigned int j) const;
00333 
00340     unsigned int row_length (const unsigned int row) const;
00341 
00361     unsigned int n_nonzero_elements () const;
00362 
00372     void print (std::ostream &out) const;
00373 
00383     void print_gnuplot (std::ostream &out) const;
00384 
00391     DeclException4 (ExcIncompatibleRowNumbers,
00392                     int, int, int, int,
00393                     << "The blocks [" << arg1 << ',' << arg2 << "] and ["
00394                     << arg3 << ',' << arg4 << "] have differing row numbers.");
00398     DeclException4 (ExcIncompatibleColNumbers,
00399                     int, int, int, int,
00400                     << "The blocks [" << arg1 << ',' << arg2 << "] and ["
00401                     << arg3 << ',' << arg4 << "] have differing column numbers.");
00405     DeclException0 (ExcInvalidConstructorCall);
00407 
00408   protected:
00409 
00413     unsigned int rows;
00414 
00418     unsigned int columns;
00419 
00423     Table<2,SmartPointer<SparsityPatternBase, BlockSparsityPatternBase<SparsityPatternBase> > > sub_objects;
00424 
00431     BlockIndices    row_indices;
00432 
00439     BlockIndices    column_indices;
00440 
00441   private:
00448     std::vector<unsigned int> counter_within_block;
00449 
00456     std::vector<std::vector<unsigned int> > block_column_indices;
00457 
00464     template <typename number>
00465     friend class BlockSparseMatrix;
00466 };
00467 
00468 
00469 
00480 class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
00481 {
00482   public:
00483 
00494     BlockSparsityPattern ();
00495 
00504     BlockSparsityPattern (const unsigned int n_rows,
00505                           const unsigned int n_columns);
00506 
00511     void reinit (const unsigned int n_block_rows,
00512                  const unsigned int n_block_columns);
00513 
00536     void reinit (const BlockIndices& row_indices,
00537                  const BlockIndices& col_indices,
00538                  const std::vector<std::vector<unsigned int> >& row_lengths);
00539 
00540 
00547     bool is_compressed () const;
00548 
00554     std::size_t memory_consumption () const;
00555 
00567     void copy_from (const BlockCompressedSparsityPattern &csp);
00568 
00580     void copy_from (const BlockCompressedSetSparsityPattern &csp);
00581 
00593     void copy_from (const BlockCompressedSimpleSparsityPattern &csp);
00594 
00595 
00596 
00597 };
00598 
00599 
00600 
00651 class BlockCompressedSparsityPattern : public BlockSparsityPatternBase<CompressedSparsityPattern>
00652 {
00653   public:
00654 
00665     BlockCompressedSparsityPattern ();
00666 
00675     BlockCompressedSparsityPattern (const unsigned int n_rows,
00676                                     const unsigned int n_columns);
00677 
00689     BlockCompressedSparsityPattern (const std::vector<unsigned int>& row_block_sizes,
00690                                     const std::vector<unsigned int>& col_block_sizes);
00691 
00698     BlockCompressedSparsityPattern (const BlockIndices& row_indices,
00699                                     const BlockIndices& col_indices);
00700 
00716     void reinit (const std::vector< unsigned int > &row_block_sizes,
00717                  const std::vector< unsigned int > &col_block_sizes);
00718 
00729     void reinit (const BlockIndices& row_indices, const BlockIndices& col_indices);
00730 
00736     using BlockSparsityPatternBase<CompressedSparsityPattern>::reinit;
00737 };
00738 
00739 
00748 typedef BlockCompressedSparsityPattern CompressedBlockSparsityPattern;
00749 
00750 
00751 
00768 class BlockCompressedSetSparsityPattern : public BlockSparsityPatternBase<CompressedSetSparsityPattern>
00769 {
00770   public:
00771 
00782     BlockCompressedSetSparsityPattern ();
00783 
00792     BlockCompressedSetSparsityPattern (const unsigned int n_rows,
00793                                        const unsigned int n_columns);
00794 
00806     BlockCompressedSetSparsityPattern (const std::vector<unsigned int>& row_block_sizes,
00807                                        const std::vector<unsigned int>& col_block_sizes);
00808 
00815     BlockCompressedSetSparsityPattern (const BlockIndices& row_indices,
00816                                        const BlockIndices& col_indices);
00817 
00833     void reinit (const std::vector< unsigned int > &row_block_sizes,
00834                  const std::vector< unsigned int > &col_block_sizes);
00835 
00846     void reinit (const BlockIndices& row_indices, const BlockIndices& col_indices);
00847 
00853     using BlockSparsityPatternBase<CompressedSetSparsityPattern>::reinit;
00854 };
00855 
00856 
00857 
00858 
00859 
00878 class BlockCompressedSimpleSparsityPattern : public BlockSparsityPatternBase<CompressedSimpleSparsityPattern>
00879 {
00880   public:
00881 
00892     BlockCompressedSimpleSparsityPattern ();
00893 
00902     BlockCompressedSimpleSparsityPattern (const unsigned int n_rows,
00903                                           const unsigned int n_columns);
00904 
00916     BlockCompressedSimpleSparsityPattern (const std::vector<unsigned int>& row_block_sizes,
00917                                     const std::vector<unsigned int>& col_block_sizes);
00918 
00932     BlockCompressedSimpleSparsityPattern (const std::vector<IndexSet> & partitioning);
00933 
00948     void reinit (const std::vector< unsigned int > &row_block_sizes,
00949                  const std::vector< unsigned int > &col_block_sizes);
00950 
00958     void reinit(const std::vector<IndexSet> & partitioning);
00959 
00965     using BlockSparsityPatternBase<CompressedSimpleSparsityPattern>::reinit;
00966 };
00967 
00968 
00969 
00970 
00971 #ifdef DEAL_II_USE_TRILINOS
00972 
00973 
00993 namespace TrilinosWrappers
00994 {
00995   class BlockSparsityPattern :
00996     public ::BlockSparsityPatternBase<SparsityPattern>
00997   {
00998     public:
00999 
01010       BlockSparsityPattern ();
01011 
01020       BlockSparsityPattern (const unsigned int n_rows,
01021                             const unsigned int n_columns);
01022 
01034       BlockSparsityPattern (const std::vector<unsigned int>& row_block_sizes,
01035                             const std::vector<unsigned int>& col_block_sizes);
01036 
01051       BlockSparsityPattern (const std::vector<Epetra_Map>& parallel_partitioning);
01052 
01064       BlockSparsityPattern (const std::vector<IndexSet>& parallel_partitioning,
01065                             const MPI_Comm &communicator = MPI_COMM_WORLD);
01066 
01082       void reinit (const std::vector< unsigned int > &row_block_sizes,
01083                    const std::vector< unsigned int > &col_block_sizes);
01084 
01092       void reinit (const std::vector<Epetra_Map>& parallel_partitioning);
01093 
01100       void reinit (const std::vector<IndexSet>& parallel_partitioning,
01101                    const MPI_Comm             & communicator = MPI_COMM_WORLD);
01102 
01103 
01109     using BlockSparsityPatternBase<SparsityPattern>::reinit;
01110   };
01111 }
01112 
01113 #endif
01114 
01115 
01117 /*---------------------- Template functions -----------------------------------*/
01118 
01119 
01120 
01121 template <class SparsityPatternBase>
01122 inline
01123 SparsityPatternBase &
01124 BlockSparsityPatternBase<SparsityPatternBase>::block (const unsigned int row,
01125                                                       const unsigned int column)
01126 {
01127   Assert (row<rows, ExcIndexRange(row,0,rows));
01128   Assert (column<columns, ExcIndexRange(column,0,columns));
01129   return *sub_objects[row][column];
01130 }
01131 
01132 
01133 
01134 template <class SparsityPatternBase>
01135 inline
01136 const SparsityPatternBase &
01137 BlockSparsityPatternBase<SparsityPatternBase>::block (const unsigned int row,
01138                                                       const unsigned int column) const
01139 {
01140   Assert (row<rows, ExcIndexRange(row,0,rows));
01141   Assert (column<columns, ExcIndexRange(column,0,columns));
01142   return *sub_objects[row][column];
01143 }
01144 
01145 
01146 
01147 template <class SparsityPatternBase>
01148 inline
01149 const BlockIndices &
01150 BlockSparsityPatternBase<SparsityPatternBase>::get_row_indices () const
01151 {
01152   return row_indices;
01153 }
01154 
01155 
01156 
01157 template <class SparsityPatternBase>
01158 inline
01159 const BlockIndices &
01160 BlockSparsityPatternBase<SparsityPatternBase>::get_column_indices () const
01161 {
01162   return column_indices;
01163 }
01164 
01165 
01166 
01167 template <class SparsityPatternBase>
01168 inline
01169 void
01170 BlockSparsityPatternBase<SparsityPatternBase>::add (const unsigned int i,
01171                                                     const unsigned int j)
01172 {
01173                                    // if you get an error here, are
01174                                    // you sure you called
01175                                    // <tt>collect_sizes()</tt> before?
01176   const std::pair<unsigned int,unsigned int>
01177     row_index = row_indices.global_to_local (i),
01178     col_index = column_indices.global_to_local (j);
01179   sub_objects[row_index.first][col_index.first]->add (row_index.second,
01180                                                       col_index.second);
01181 }
01182 
01183 
01184 
01185 template <class SparsityPatternBase>
01186 template <typename ForwardIterator>
01187 void
01188 BlockSparsityPatternBase<SparsityPatternBase>::add_entries (const unsigned int row,
01189                                                             ForwardIterator    begin,
01190                                                             ForwardIterator    end,
01191                                                             const bool         indices_are_sorted)
01192 {
01193                                    // Resize scratch arrays
01194   if (block_column_indices.size() < this->n_block_cols())
01195     {
01196       block_column_indices.resize (this->n_block_cols());
01197       counter_within_block.resize (this->n_block_cols());
01198     }
01199 
01200   const unsigned int n_cols = static_cast<unsigned int>(end - begin);
01201 
01202                                    // Resize sub-arrays to n_cols. This
01203                                    // is a bit wasteful, but we resize
01204                                    // only a few times (then the maximum
01205                                    // row length won't increase that
01206                                    // much any more). At least we know
01207                                    // that all arrays are going to be of
01208                                    // the same size, so we can check
01209                                    // whether the size of one is large
01210                                    // enough before actually going
01211                                    // through all of them.
01212   if (block_column_indices[0].size() < n_cols)
01213     for (unsigned int i=0; i<this->n_block_cols(); ++i)
01214       block_column_indices[i].resize(n_cols);
01215 
01216                                    // Reset the number of added elements
01217                                    // in each block to zero.
01218   for (unsigned int i=0; i<this->n_block_cols(); ++i)
01219     counter_within_block[i] = 0;
01220 
01221                                    // Go through the column indices to
01222                                    // find out which portions of the
01223                                    // values should be set in which
01224                                    // block of the matrix. We need to
01225                                    // touch all the data, since we can't
01226                                    // be sure that the data of one block
01227                                    // is stored contiguously (in fact,
01228                                    // indices will be intermixed when it
01229                                    // comes from an element matrix).
01230   for (ForwardIterator it = begin; it != end; ++it)
01231     {
01232       const unsigned int col = *it;
01233 
01234       const std::pair<unsigned int, unsigned int>
01235         col_index = this->column_indices.global_to_local(col);
01236 
01237       const unsigned int local_index = counter_within_block[col_index.first]++;
01238 
01239       block_column_indices[col_index.first][local_index] = col_index.second;
01240     }
01241 
01242 #ifdef DEBUG
01243                                    // If in debug mode, do a check whether
01244                                    // the right length has been obtained.
01245   unsigned int length = 0;
01246   for (unsigned int i=0; i<this->n_block_cols(); ++i)
01247     length += counter_within_block[i];
01248   Assert (length == n_cols, ExcInternalError());
01249 #endif
01250 
01251                                    // Now we found out about where the
01252                                    // individual columns should start and
01253                                    // where we should start reading out
01254                                    // data. Now let's write the data into
01255                                    // the individual blocks!
01256   const std::pair<unsigned int,unsigned int>
01257     row_index = this->row_indices.global_to_local (row);
01258   for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
01259     {
01260       if (counter_within_block[block_col] == 0)
01261         continue;
01262       sub_objects[row_index.first][block_col]->
01263         add_entries (row_index.second,
01264                      block_column_indices[block_col].begin(),
01265                      block_column_indices[block_col].begin()+counter_within_block[block_col],
01266                      indices_are_sorted);
01267     }
01268 }
01269 
01270 
01271 
01272 template <class SparsityPatternBase>
01273 inline
01274 bool
01275 BlockSparsityPatternBase<SparsityPatternBase>::exists (const unsigned int i,
01276                                                        const unsigned int j) const
01277 {
01278                                    // if you get an error here, are
01279                                    // you sure you called
01280                                    // <tt>collect_sizes()</tt> before?
01281   const std::pair<unsigned int,unsigned int>
01282     row_index = row_indices.global_to_local (i),
01283     col_index = column_indices.global_to_local (j);
01284   return sub_objects[row_index.first][col_index.first]->exists (row_index.second,
01285                                                                 col_index.second);
01286 }
01287 
01288 
01289 
01290 template <class SparsityPatternBase>
01291 inline
01292 unsigned int
01293 BlockSparsityPatternBase<SparsityPatternBase>::
01294 row_length (const unsigned int row) const
01295 {
01296   const std::pair<unsigned int,unsigned int>
01297     row_index = row_indices.global_to_local (row);
01298 
01299   unsigned int c = 0;
01300 
01301   for (unsigned int b=0; b<rows; ++b)
01302     c += sub_objects[row_index.first][b]->row_length (row_index.second);
01303 
01304   return c;
01305 }
01306 
01307 
01308 
01309 template <class SparsityPatternBase>
01310 inline
01311 unsigned int
01312 BlockSparsityPatternBase<SparsityPatternBase>::n_block_cols () const
01313 {
01314   return columns;
01315 }
01316 
01317 
01318 
01319 template <class SparsityPatternBase>
01320 inline
01321 unsigned int
01322 BlockSparsityPatternBase<SparsityPatternBase>::n_block_rows () const
01323 {
01324   return rows;
01325 }
01326 
01327 
01328 inline
01329 void
01330 BlockSparsityPattern::reinit (
01331   const unsigned int n_block_rows,
01332   const unsigned int n_block_columns)
01333 {
01334   BlockSparsityPatternBase<SparsityPattern>::reinit (
01335     n_block_rows, n_block_columns);
01336 }
01337 
01338 
01339 DEAL_II_NAMESPACE_CLOSE
01340 
01341 #endif
01342 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Thu May 17 2012 20:04:40 by doxygen 1.7.3