Reference documentation for deal.II version 8.1.0
block_sparsity_pattern.h
1 // ---------------------------------------------------------------------
2 // @f$Id: block_sparsity_pattern.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2000 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__block_sparsity_pattern_h
18 #define __deal2__block_sparsity_pattern_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/table.h>
24 #include <deal.II/base/subscriptor.h>
25 #include <deal.II/base/smartpointer.h>
26 #include <deal.II/lac/sparsity_pattern.h>
27 #include <deal.II/lac/trilinos_sparsity_pattern.h>
28 #include <deal.II/lac/compressed_sparsity_pattern.h>
29 #include <deal.II/lac/compressed_set_sparsity_pattern.h>
30 #include <deal.II/lac/compressed_simple_sparsity_pattern.h>
31 #include <deal.II/lac/block_indices.h>
32 
34 
35 
36 template <typename number> class BlockSparseMatrix;
41 #ifdef DEAL_II_WITH_TRILINOS
42 namespace TrilinosWrappers
43 {
45 }
46 #endif
47 
86 template <class SparsityPatternBase>
88 {
89 public:
94 
108 
120 
129  BlockSparsityPatternBase (const size_type n_block_rows,
130  const size_type n_block_columns);
131 
144 
149 
176  void reinit (const size_type n_block_rows,
177  const size_type n_block_columns);
178 
188 
201  void collect_sizes ();
202 
207  SparsityPatternBase &
208  block (const size_type row,
209  const size_type column);
210 
211 
217  const SparsityPatternBase &
218  block (const size_type row,
219  const size_type column) const;
220 
227  const BlockIndices &
228  get_row_indices () const;
229 
236  const BlockIndices &
237  get_column_indices () const;
238 
246  void compress ();
247 
252  size_type n_block_rows () const;
253 
258  size_type n_block_cols () const;
259 
271  bool empty () const;
272 
281  size_type max_entries_per_row () const;
282 
295  void add (const size_type i, const size_type j);
296 
312  template <typename ForwardIterator>
313  void add_entries (const size_type row,
314  ForwardIterator begin,
315  ForwardIterator end,
316  const bool indices_are_sorted = false);
317 
326  size_type n_rows () const;
327 
336  size_type n_cols () const;
337 
342  bool exists (const size_type i, const size_type j) const;
343 
350  unsigned int row_length (const size_type row) const;
351 
371  size_type n_nonzero_elements () const;
372 
382  void print (std::ostream &out) const;
383 
393  void print_gnuplot (std::ostream &out) const;
394 
401  DeclException4 (ExcIncompatibleRowNumbers,
402  int, int, int, int,
403  << "The blocks [" << arg1 << ',' << arg2 << "] and ["
404  << arg3 << ',' << arg4 << "] have differing row numbers.");
408  DeclException4 (ExcIncompatibleColNumbers,
409  int, int, int, int,
410  << "The blocks [" << arg1 << ',' << arg2 << "] and ["
411  << arg3 << ',' << arg4 << "] have differing column numbers.");
417 
418 protected:
419 
423  size_type rows;
424 
428  size_type columns;
429 
434 
442 
450 
451 private:
458  std::vector<size_type > counter_within_block;
459 
466  std::vector<std::vector<size_type > > block_column_indices;
467 
474  template <typename number>
475  friend class BlockSparseMatrix;
476 };
477 
478 
479 
490 class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
491 {
492 public:
493 
505 
515  const size_type n_columns);
516 
521  void reinit (const size_type n_block_rows,
522  const size_type n_block_columns);
523 
546  void reinit (const BlockIndices &row_indices,
547  const BlockIndices &col_indices,
548  const std::vector<std::vector<unsigned int> > &row_lengths);
549 
550 
557  bool is_compressed () const;
558 
564  std::size_t memory_consumption () const;
565 
577  void copy_from (const BlockCompressedSparsityPattern &csp);
578 
591 
604 
605 
606 
607 };
608 
609 
610 
661 class BlockCompressedSparsityPattern : public BlockSparsityPatternBase<CompressedSparsityPattern>
662 {
663 public:
664 
676 
686  const size_type n_columns);
687 
699  BlockCompressedSparsityPattern (const std::vector<size_type> &row_block_sizes,
700  const std::vector<size_type> &col_block_sizes);
701 
709  const BlockIndices &col_indices);
710 
726  void reinit (const std::vector<size_type> &row_block_sizes,
727  const std::vector<size_type> &col_block_sizes);
728 
739  void reinit (const BlockIndices &row_indices, const BlockIndices &col_indices);
740 
747 };
748 
749 
758 typedef BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED;
759 
760 
761 
778 class BlockCompressedSetSparsityPattern : public BlockSparsityPatternBase<CompressedSetSparsityPattern>
779 {
780 public:
781 
793 
803  const size_type n_columns);
804 
816  BlockCompressedSetSparsityPattern (const std::vector<size_type> &row_block_sizes,
817  const std::vector<size_type> &col_block_sizes);
818 
826  const BlockIndices &col_indices);
827 
843  void reinit (const std::vector<size_type> &row_block_sizes,
844  const std::vector<size_type> &col_block_sizes);
845 
856  void reinit (const BlockIndices &row_indices, const BlockIndices &col_indices);
857 
864 };
865 
866 
867 
868 
869 
888 class BlockCompressedSimpleSparsityPattern : public BlockSparsityPatternBase<CompressedSimpleSparsityPattern>
889 {
890 public:
891 
903 
913  const size_type n_columns);
914 
926  BlockCompressedSimpleSparsityPattern (const std::vector<size_type> &row_block_sizes,
927  const std::vector<size_type> &col_block_sizes);
928 
942  BlockCompressedSimpleSparsityPattern (const std::vector<IndexSet> &partitioning);
943 
958  void reinit (const std::vector<size_type> &row_block_sizes,
959  const std::vector<size_type> &col_block_sizes);
960 
968  void reinit(const std::vector<IndexSet> &partitioning);
969 
975  size_type column_number (const size_type row,
976  const unsigned int index) const;
977 
984 };
985 
986 
987 
988 
989 #ifdef DEAL_II_WITH_TRILINOS
990 
991 
1011 namespace TrilinosWrappers
1012 {
1014  public ::BlockSparsityPatternBase<SparsityPattern>
1015  {
1016  public:
1017 
1029 
1039  const size_type n_columns);
1040 
1052  BlockSparsityPattern (const std::vector<size_type> &row_block_sizes,
1053  const std::vector<size_type> &col_block_sizes);
1054 
1069  BlockSparsityPattern (const std::vector<Epetra_Map> &parallel_partitioning);
1070 
1082  BlockSparsityPattern (const std::vector<IndexSet> &parallel_partitioning,
1083  const MPI_Comm &communicator = MPI_COMM_WORLD);
1084 
1100  void reinit (const std::vector<size_type> &row_block_sizes,
1101  const std::vector<size_type> &col_block_sizes);
1102 
1110  void reinit (const std::vector<Epetra_Map> &parallel_partitioning);
1111 
1118  void reinit (const std::vector<IndexSet> &parallel_partitioning,
1119  const MPI_Comm &communicator = MPI_COMM_WORLD);
1120 
1121 
1128  };
1129 }
1130 
1131 #endif
1132 
1133 
1135 /*---------------------- Template functions -----------------------------------*/
1136 
1137 
1138 
1139 template <class SparsityPatternBase>
1140 inline
1141 SparsityPatternBase &
1143  const size_type column)
1144 {
1145  Assert (row<rows, ExcIndexRange(row,0,rows));
1146  Assert (column<columns, ExcIndexRange(column,0,columns));
1147  return *sub_objects[row][column];
1148 }
1149 
1150 
1151 
1152 template <class SparsityPatternBase>
1153 inline
1154 const SparsityPatternBase &
1156  const size_type column) const
1157 {
1158  Assert (row<rows, ExcIndexRange(row,0,rows));
1159  Assert (column<columns, ExcIndexRange(column,0,columns));
1160  return *sub_objects[row][column];
1161 }
1162 
1163 
1164 
1165 template <class SparsityPatternBase>
1166 inline
1167 const BlockIndices &
1169 {
1170  return row_indices;
1171 }
1172 
1173 
1174 
1175 template <class SparsityPatternBase>
1176 inline
1177 const BlockIndices &
1179 {
1180  return column_indices;
1181 }
1182 
1183 
1184 
1185 template <class SparsityPatternBase>
1186 inline
1187 void
1189  const size_type j)
1190 {
1191  // if you get an error here, are
1192  // you sure you called
1193  // <tt>collect_sizes()</tt> before?
1194  const std::pair<size_type,size_type>
1195  row_index = row_indices.global_to_local (i),
1196  col_index = column_indices.global_to_local (j);
1197  sub_objects[row_index.first][col_index.first]->add (row_index.second,
1198  col_index.second);
1199 }
1200 
1201 
1202 
1203 template <class SparsityPatternBase>
1204 template <typename ForwardIterator>
1205 void
1207  ForwardIterator begin,
1208  ForwardIterator end,
1209  const bool indices_are_sorted)
1210 {
1211  // Resize scratch arrays
1212  if (block_column_indices.size() < this->n_block_cols())
1213  {
1214  block_column_indices.resize (this->n_block_cols());
1215  counter_within_block.resize (this->n_block_cols());
1216  }
1217 
1218  const size_type n_cols = static_cast<size_type>(end - begin);
1219 
1220  // Resize sub-arrays to n_cols. This
1221  // is a bit wasteful, but we resize
1222  // only a few times (then the maximum
1223  // row length won't increase that
1224  // much any more). At least we know
1225  // that all arrays are going to be of
1226  // the same size, so we can check
1227  // whether the size of one is large
1228  // enough before actually going
1229  // through all of them.
1230  if (block_column_indices[0].size() < n_cols)
1231  for (size_type i=0; i<this->n_block_cols(); ++i)
1232  block_column_indices[i].resize(n_cols);
1233 
1234  // Reset the number of added elements
1235  // in each block to zero.
1236  for (size_type i=0; i<this->n_block_cols(); ++i)
1237  counter_within_block[i] = 0;
1238 
1239  // Go through the column indices to
1240  // find out which portions of the
1241  // values should be set in which
1242  // block of the matrix. We need to
1243  // touch all the data, since we can't
1244  // be sure that the data of one block
1245  // is stored contiguously (in fact,
1246  // indices will be intermixed when it
1247  // comes from an element matrix).
1248  for (ForwardIterator it = begin; it != end; ++it)
1249  {
1250  const size_type col = *it;
1251 
1252  const std::pair<size_type , size_type>
1253  col_index = this->column_indices.global_to_local(col);
1254 
1255  const size_type local_index = counter_within_block[col_index.first]++;
1256 
1257  block_column_indices[col_index.first][local_index] = col_index.second;
1258  }
1259 
1260 #ifdef DEBUG
1261  // If in debug mode, do a check whether
1262  // the right length has been obtained.
1263  size_type length = 0;
1264  for (size_type i=0; i<this->n_block_cols(); ++i)
1265  length += counter_within_block[i];
1266  Assert (length == n_cols, ExcInternalError());
1267 #endif
1268 
1269  // Now we found out about where the
1270  // individual columns should start and
1271  // where we should start reading out
1272  // data. Now let's write the data into
1273  // the individual blocks!
1274  const std::pair<size_type , size_type>
1275  row_index = this->row_indices.global_to_local (row);
1276  for (size_type block_col=0; block_col<n_block_cols(); ++block_col)
1277  {
1278  if (counter_within_block[block_col] == 0)
1279  continue;
1280  sub_objects[row_index.first][block_col]->
1281  add_entries (row_index.second,
1282  block_column_indices[block_col].begin(),
1283  block_column_indices[block_col].begin()+counter_within_block[block_col],
1284  indices_are_sorted);
1285  }
1286 }
1287 
1288 
1289 
1290 template <class SparsityPatternBase>
1291 inline
1292 bool
1294  const size_type j) const
1295 {
1296  // if you get an error here, are
1297  // you sure you called
1298  // <tt>collect_sizes()</tt> before?
1299  const std::pair<size_type , size_type>
1300  row_index = row_indices.global_to_local (i),
1301  col_index = column_indices.global_to_local (j);
1302  return sub_objects[row_index.first][col_index.first]->exists (row_index.second,
1303  col_index.second);
1304 }
1305 
1306 
1307 
1308 template <class SparsityPatternBase>
1309 inline
1310 unsigned int
1312 row_length (const size_type row) const
1313 {
1314  const std::pair<size_type , size_type>
1315  row_index = row_indices.global_to_local (row);
1316 
1317  unsigned int c = 0;
1318 
1319  for (size_type b=0; b<rows; ++b)
1320  c += sub_objects[row_index.first][b]->row_length (row_index.second);
1321 
1322  return c;
1323 }
1324 
1325 
1326 
1327 template <class SparsityPatternBase>
1328 inline
1331 {
1332  return columns;
1333 }
1334 
1335 
1336 
1337 template <class SparsityPatternBase>
1338 inline
1341 {
1342  return rows;
1343 }
1344 
1345 
1346 inline
1349  const unsigned int index) const
1350 {
1351  // .first= ith block, .second = jth row in that block
1352  const std::pair<size_type ,size_type >
1353  row_index = row_indices.global_to_local (row);
1354 
1355  Assert(index<row_length(row), ExcIndexRange(index, 0, row_length(row)));
1356 
1357  size_type c = 0;
1358  size_type block_columns = 0; //sum of n_cols for all blocks to the left
1359  for (unsigned int b=0; b<columns; ++b)
1360  {
1361  unsigned int rowlen = sub_objects[row_index.first][b]->row_length (row_index.second);
1362  if (index<c+rowlen)
1363  return block_columns+sub_objects[row_index.first][b]->column_number(row_index.second, index-c);
1364  c += rowlen;
1365  block_columns += sub_objects[row_index.first][b]->n_cols();
1366  }
1367 
1368  Assert(false, ExcInternalError());
1369  return 0;
1370 }
1371 
1372 
1373 inline
1374 void
1376  const size_type n_block_rows,
1377  const size_type n_block_columns)
1378 {
1380  n_block_rows, n_block_columns);
1381 }
1382 
1383 
1384 DEAL_II_NAMESPACE_CLOSE
1385 
1386 #endif
size_type column_number(const size_type row, const unsigned int index) const
void print(std::ostream &out) const
static const size_type invalid_entry
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:55
unsigned int row_length(const size_type row) const
DeclException4(ExcIncompatibleRowNumbers, int, int, int, int,<< "The blocks ["<< arg1<< ','<< arg2<< "] and ["<< arg3<< ','<< arg4<< "] have differing row numbers.")
void add(const size_type i, const size_type j)
std::size_t memory_consumption() const
bool is_compressed() const
size_type n_cols() const
SparsityPatternBase & block(const size_type row, const size_type column)
void reinit(const size_type n_block_rows, const size_type n_block_columns)
void reinit(const size_type n_block_rows, const size_type n_block_columns)
Table< 2, SmartPointer< SparsityPatternBase, BlockSparsityPatternBase< SparsityPatternBase > > > sub_objects
unsigned int global_dof_index
Definition: types.h:100
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
size_type max_entries_per_row() const
#define Assert(cond, exc)
Definition: exceptions.h:299
types::global_dof_index size_type
const BlockIndices & get_column_indices() const
BlockSparsityPatternBase & operator=(const BlockSparsityPatternBase &)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
::ExceptionBase & ExcInvalidConstructorCall()
std::vector< size_type > counter_within_block
std::vector< std::vector< size_type > > block_column_indices
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
const BlockIndices & get_row_indices() const
size_type n_rows() const
size_type n_nonzero_elements() const
void print_gnuplot(std::ostream &out) const
DeclException0(ExcInvalidConstructorCall)
void copy_from(const BlockCompressedSparsityPattern &csp)
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
::ExceptionBase & ExcInternalError()
bool exists(const size_type i, const size_type j) const
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
static const size_type invalid_entry