Reference documentation for deal.II version 8.4.1
block_sparsity_pattern.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2015 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/vector_slice.h>
18 #include <deal.II/base/memory_consumption.h>
19 #include <deal.II/lac/block_sparsity_pattern.h>
20 
21 DEAL_II_NAMESPACE_OPEN
22 
23 
24 template <class SparsityPatternBase>
26  :
27  rows (0),
28  columns (0)
29 {}
30 
31 
32 
33 template <class SparsityPatternBase>
35 BlockSparsityPatternBase (const size_type n_block_rows,
36  const size_type n_block_columns)
37  :
38  rows (0),
39  columns (0)
40 {
41  reinit (n_block_rows, n_block_columns);
42 }
43 
44 
45 
46 template <class SparsityPatternBase>
49  :
50  Subscriptor ()
51 {
52  (void)s;
53  Assert(s.rows==0, ExcInvalidConstructorCall());
54  Assert(s.columns==0, ExcInvalidConstructorCall());
55 
56  rows = 0;
57  columns=0;
58 }
59 
60 
61 
62 template <class SparsityPatternBase>
64 {
65  // clear all memory
66  reinit (0,0);
67 }
68 
69 
70 
71 template <class SparsityPatternBase>
72 void
74 reinit (const size_type n_block_rows,
75  const size_type n_block_columns)
76 {
77  // delete previous content and
78  // clean the sub_objects array
79  // completely
80  for (size_type i=0; i<rows; ++i)
81  for (size_type j=0; j<columns; ++j)
82  {
83  SparsityPatternBase *sp = sub_objects[i][j];
84  sub_objects[i][j] = 0;
85  delete sp;
86  };
87  sub_objects.reinit (0,0);
88 
89  // then set new sizes
90  rows = n_block_rows;
91  columns = n_block_columns;
92  sub_objects.reinit (rows, columns);
93 
94  // allocate new objects
95  for (size_type i=0; i<rows; ++i)
96  for (size_type j=0; j<columns; ++j)
97  {
98  SparsityPatternBase *p = new SparsityPatternBase;
99  sub_objects[i][j] = p;
100  }
101 }
102 
103 
104 template <class SparsityPatternBase>
108 {
109  Assert (rows == bsp.rows, ExcDimensionMismatch(rows, bsp.rows));
110  Assert (columns == bsp.columns, ExcDimensionMismatch(columns, bsp.columns));
111  // copy objects
112  for (size_type i=0; i<rows; ++i)
113  for (size_type j=0; j<columns; ++j)
114  *sub_objects[i][j] = *bsp.sub_objects[i][j];
115  // update index objects
116  collect_sizes ();
117 
118  return *this;
119 }
120 
121 
122 
123 template <class SparsityPatternBase>
124 void
126 {
127  std::vector<size_type > row_sizes (rows);
128  std::vector<size_type > col_sizes (columns);
129 
130  // first find out the row sizes
131  // from the first block column
132  for (size_type r=0; r<rows; ++r)
133  row_sizes[r] = sub_objects[r][0]->n_rows();
134  // then check that the following
135  // block columns have the same
136  // sizes
137  for (size_type c=1; c<columns; ++c)
138  for (size_type r=0; r<rows; ++r)
139  Assert (row_sizes[r] == sub_objects[r][c]->n_rows(),
140  ExcIncompatibleRowNumbers (r,0,r,c));
141 
142  // finally initialize the row
143  // indices with this array
144  row_indices.reinit (row_sizes);
145 
146 
147  // then do the same with the columns
148  for (size_type c=0; c<columns; ++c)
149  col_sizes[c] = sub_objects[0][c]->n_cols();
150  for (size_type r=1; r<rows; ++r)
151  for (size_type c=0; c<columns; ++c)
152  Assert (col_sizes[c] == sub_objects[r][c]->n_cols(),
153  ExcIncompatibleRowNumbers (0,c,r,c));
154 
155  // finally initialize the row
156  // indices with this array
157  column_indices.reinit (col_sizes);
158 }
159 
160 
161 
162 template <class SparsityPatternBase>
163 void
165 {
166  for (size_type i=0; i<rows; ++i)
167  for (size_type j=0; j<columns; ++j)
168  sub_objects[i][j]->compress ();
169 }
170 
171 
172 
173 template <class SparsityPatternBase>
174 bool
176 {
177  for (size_type i=0; i<rows; ++i)
178  for (size_type j=0; j<columns; ++j)
179  if (sub_objects[i][j]->empty () == false)
180  return false;
181  return true;
182 }
183 
184 
185 
186 template <class SparsityPatternBase>
189 {
190  size_type max_entries = 0;
191  for (size_type block_row=0; block_row<rows; ++block_row)
192  {
193  size_type this_row = 0;
194  for (size_type c=0; c<columns; ++c)
195  this_row += sub_objects[block_row][c]->max_entries_per_row ();
196 
197  if (this_row > max_entries)
198  max_entries = this_row;
199  };
200  return max_entries;
201 }
202 
203 
204 
205 template <class SparsityPatternBase>
208 {
209  // only count in first column, since
210  // all rows should be equivalent
211  size_type count = 0;
212  for (size_type r=0; r<rows; ++r)
213  count += sub_objects[r][0]->n_rows();
214  return count;
215 }
216 
217 
218 
219 template <class SparsityPatternBase>
222 {
223  // only count in first row, since
224  // all rows should be equivalent
225  size_type count = 0;
226  for (size_type c=0; c<columns; ++c)
227  count += sub_objects[0][c]->n_cols();
228  return count;
229 }
230 
231 
232 
233 template <class SparsityPatternBase>
236 {
237  size_type count = 0;
238  for (size_type i=0; i<rows; ++i)
239  for (size_type j=0; j<columns; ++j)
240  count += sub_objects[i][j]->n_nonzero_elements ();
241  return count;
242 }
243 
244 
245 
246 template <class SparsityPatternBase>
247 void
249 {
250  size_type k=0;
251  for (size_type ib=0; ib<n_block_rows(); ++ib)
252  {
253  for (size_type i=0; i<block(ib,0).n_rows(); ++i)
254  {
255  out << '[' << i+k;
256  size_type l=0;
257  for (size_type jb=0; jb<n_block_cols(); ++jb)
258  {
259  const SparsityPatternBase &b = block(ib,jb);
260  for (size_type j=0; j<b.n_cols(); ++j)
261  if (b.exists(i,j))
262  out << ',' << l+j;
263  l += b.n_cols();
264  }
265  out << ']' << std::endl;
266  }
267  k += block(ib,0).n_rows();
268  }
269 }
270 
271 
272 template <>
273 void
275 {
276  size_type k=0;
277  for (size_type ib=0; ib<n_block_rows(); ++ib)
278  {
279  for (size_type i=0; i<block(ib,0).n_rows(); ++i)
280  {
281  out << '[' << i+k;
282  size_type l=0;
283  for (size_type jb=0; jb<n_block_cols(); ++jb)
284  {
285  const DynamicSparsityPattern &b = block(ib,jb);
286  if (b.row_index_set().size()==0 || b.row_index_set().is_element(i))
287  for (size_type j=0; j<b.n_cols(); ++j)
288  if (b.exists(i,j))
289  out << ',' << l+j;
290  l += b.n_cols();
291  }
292  out << ']' << std::endl;
293  }
294  k += block(ib,0).n_rows();
295  }
296 }
297 
298 
299 template <class SparsityPatternBase>
300 void
302 {
303  size_type k=0;
304  for (size_type ib=0; ib<n_block_rows(); ++ib)
305  {
306  for (size_type i=0; i<block(ib,0).n_rows(); ++i)
307  {
308  size_type l=0;
309  for (size_type jb=0; jb<n_block_cols(); ++jb)
310  {
311  const SparsityPatternBase &b = block(ib,jb);
312  for (size_type j=0; j<b.n_cols(); ++j)
313  if (b.exists(i,j))
314  out << l+j << " " << -static_cast<signed int>(i+k) << std::endl;
315  l += b.n_cols();
316  }
317  }
318  k += block(ib,0).n_rows();
319  }
320 }
321 
322 
323 
325 {}
326 
327 
328 
330  const size_type n_columns)
331  :
333  n_columns)
334 {}
335 
336 
337 void
339  const BlockIndices &rows,
340  const BlockIndices &cols,
341  const std::vector<std::vector<unsigned int> > &row_lengths)
342 {
343  AssertDimension (row_lengths.size(), cols.size());
344 
345  this->reinit(rows.size(), cols.size());
346  for (size_type j=0; j<cols.size(); ++j)
347  for (size_type i=0; i<rows.size(); ++i)
348  {
349  const size_type start = rows.local_to_global(i, 0);
350  const size_type length = rows.block_size(i);
351 
352  if (row_lengths[j].size()==1)
353  block(i,j).reinit(rows.block_size(i),
354  cols.block_size(j), row_lengths[j][0]);
355  else
356  {
358  block_rows(row_lengths[j], start, length);
359  block(i,j).reinit(rows.block_size(i),
360  cols.block_size(j),
361  block_rows);
362  }
363  }
364  this->collect_sizes();
365  Assert (this->row_indices == rows, ExcInternalError());
366  Assert (this->column_indices == cols, ExcInternalError());
367 }
368 
369 
370 bool
372 {
373  for (size_type i=0; i<rows; ++i)
374  for (size_type j=0; j<columns; ++j)
375  if (sub_objects[i][j]->is_compressed () == false)
376  return false;
377  return true;
378 }
379 
380 
381 std::size_t
383 {
384  std::size_t mem = 0;
390  for (size_type r=0; r<rows; ++r)
391  for (size_type c=0; c<columns; ++c)
393 
394  return mem;
395 }
396 
397 
398 
399 void
401 {
402  // delete old content, set block
403  // sizes anew
404  reinit (dsp.n_block_rows(), dsp.n_block_cols());
405 
406  // copy over blocks
407  for (size_type i=0; i<n_block_rows(); ++i)
408  for (size_type j=0; j<n_block_cols(); ++j)
409  block(i,j).copy_from (dsp.block(i,j));
410 
411  // and finally enquire their new
412  // sizes
413  collect_sizes();
414 }
415 
416 
417 
418 
420 {}
421 
422 
423 
426  const size_type n_columns)
427  :
429  n_columns)
430 {}
431 
432 
433 
435 BlockDynamicSparsityPattern (const std::vector<size_type> &row_indices,
436  const std::vector<size_type> &col_indices)
437  :
439  col_indices.size())
440 {
441  for (size_type i=0; i<row_indices.size(); ++i)
442  for (size_type j=0; j<col_indices.size(); ++j)
443  this->block(i,j).reinit(row_indices[i],col_indices[j]);
444  this->collect_sizes();
445 }
446 
447 
449 BlockDynamicSparsityPattern (const std::vector<IndexSet> &partitioning)
450  :
452  partitioning.size())
453 {
454  for (size_type i=0; i<partitioning.size(); ++i)
455  for (size_type j=0; j<partitioning.size(); ++j)
456  this->block(i,j).reinit(partitioning[i].size(),
457  partitioning[j].size(),
458  partitioning[i]);
459  this->collect_sizes();
460 }
461 
462 
465  const BlockIndices &col_indices)
466 {
467  reinit(row_indices, col_indices);
468 }
469 
470 
471 void
473  const std::vector<size_type> &row_block_sizes,
474  const std::vector<size_type> &col_block_sizes)
475 {
477  reinit(row_block_sizes.size(), col_block_sizes.size());
478  for (size_type i=0; i<row_block_sizes.size(); ++i)
479  for (size_type j=0; j<col_block_sizes.size(); ++j)
480  this->block(i,j).reinit(row_block_sizes[i],col_block_sizes[j]);
481  this->collect_sizes();
482 }
483 
484 void
486  const std::vector< IndexSet > &partitioning)
487 {
489  reinit(partitioning.size(), partitioning.size());
490  for (size_type i=0; i<partitioning.size(); ++i)
491  for (size_type j=0; j<partitioning.size(); ++j)
492  this->block(i,j).reinit(partitioning[i].size(),
493  partitioning[j].size(),
494  partitioning[i]);
495  this->collect_sizes();
496 }
497 
498 void
500  const BlockIndices &row_indices,
501  const BlockIndices &col_indices)
502 {
504  col_indices.size());
505  for (size_type i=0; i<row_indices.size(); ++i)
506  for (size_type j=0; j<col_indices.size(); ++j)
507  this->block(i,j).reinit(row_indices.block_size(i),
508  col_indices.block_size(j));
509  this->collect_sizes();
510 }
511 
512 
513 #ifdef DEAL_II_WITH_TRILINOS
514 namespace TrilinosWrappers
515 {
516 
518  {}
519 
520 
521 
523  BlockSparsityPattern (const size_type n_rows,
524  const size_type n_columns)
525  :
527  n_columns)
528  {}
529 
530 
531 
533  BlockSparsityPattern (const std::vector<size_type> &row_indices,
534  const std::vector<size_type> &col_indices)
535  :
536  BlockSparsityPatternBase<SparsityPattern>(row_indices.size(),
537  col_indices.size())
538  {
539  for (size_type i=0; i<row_indices.size(); ++i)
540  for (size_type j=0; j<col_indices.size(); ++j)
541  this->block(i,j).reinit(row_indices[i],col_indices[j]);
542  this->collect_sizes();
543  }
544 
545 
546 
548  BlockSparsityPattern (const std::vector<Epetra_Map> &parallel_partitioning)
549  :
551  (parallel_partitioning.size(),
552  parallel_partitioning.size())
553  {
554  for (size_type i=0; i<parallel_partitioning.size(); ++i)
555  for (size_type j=0; j<parallel_partitioning.size(); ++j)
556  this->block(i,j).reinit(parallel_partitioning[i],
557  parallel_partitioning[j]);
558  this->collect_sizes();
559  }
560 
561 
562 
564  BlockSparsityPattern (const std::vector<IndexSet> &parallel_partitioning,
565  const MPI_Comm &communicator)
566  :
568  (parallel_partitioning.size(),
569  parallel_partitioning.size())
570  {
571  for (size_type i=0; i<parallel_partitioning.size(); ++i)
572  for (size_type j=0; j<parallel_partitioning.size(); ++j)
573  this->block(i,j).reinit(parallel_partitioning[i],
574  parallel_partitioning[j],
575  communicator);
576  this->collect_sizes();
577  }
578 
579 
580 
582  BlockSparsityPattern (const std::vector<IndexSet> &row_parallel_partitioning,
583  const std::vector<IndexSet> &col_parallel_partitioning,
584  const std::vector<IndexSet> &writable_rows,
585  const MPI_Comm &communicator)
586  :
588  (row_parallel_partitioning.size(),
589  col_parallel_partitioning.size())
590  {
591  for (size_type i=0; i<row_parallel_partitioning.size(); ++i)
592  for (size_type j=0; j<col_parallel_partitioning.size(); ++j)
593  this->block(i,j).reinit(row_parallel_partitioning[i],
594  col_parallel_partitioning[j],
595  writable_rows[i],
596  communicator);
597  this->collect_sizes();
598  }
599 
600 
601 
602  void
603  BlockSparsityPattern::reinit (const std::vector<size_type> &row_block_sizes,
604  const std::vector<size_type> &col_block_sizes)
605  {
607  reinit(row_block_sizes.size(), col_block_sizes.size());
608  for (size_type i=0; i<row_block_sizes.size(); ++i)
609  for (size_type j=0; j<col_block_sizes.size(); ++j)
610  this->block(i,j).reinit(row_block_sizes[i],col_block_sizes[j]);
611  this->collect_sizes();
612  }
613 
614 
615 
616  void
617  BlockSparsityPattern::reinit (const std::vector<Epetra_Map> &parallel_partitioning)
618  {
620  reinit(parallel_partitioning.size(),
621  parallel_partitioning.size());
622  for (size_type i=0; i<parallel_partitioning.size(); ++i)
623  for (size_type j=0; j<parallel_partitioning.size(); ++j)
624  this->block(i,j).reinit(parallel_partitioning[i],
625  parallel_partitioning[j]);
626  this->collect_sizes();
627  }
628 
629 
630 
631  void
632  BlockSparsityPattern::reinit (const std::vector<IndexSet> &parallel_partitioning,
633  const MPI_Comm &communicator)
634  {
636  reinit(parallel_partitioning.size(),
637  parallel_partitioning.size());
638  for (size_type i=0; i<parallel_partitioning.size(); ++i)
639  for (size_type j=0; j<parallel_partitioning.size(); ++j)
640  this->block(i,j).reinit(parallel_partitioning[i],
641  parallel_partitioning[j],
642  communicator);
643  this->collect_sizes();
644  }
645 
646 
647 
648  void
649  BlockSparsityPattern::reinit (const std::vector<IndexSet> &row_parallel_partitioning,
650  const std::vector<IndexSet> &col_parallel_partitioning,
651  const MPI_Comm &communicator)
652  {
654  reinit(row_parallel_partitioning.size(),
655  col_parallel_partitioning.size());
656  for (size_type i=0; i<row_parallel_partitioning.size(); ++i)
657  for (size_type j=0; j<col_parallel_partitioning.size(); ++j)
658  this->block(i,j).reinit(row_parallel_partitioning[i],
659  col_parallel_partitioning[j],
660  communicator);
661  this->collect_sizes();
662  }
663 
664 
665 
666  void
667  BlockSparsityPattern::reinit (const std::vector<IndexSet> &row_parallel_partitioning,
668  const std::vector<IndexSet> &col_parallel_partitioning,
669  const std::vector<IndexSet> &writable_rows,
670  const MPI_Comm &communicator)
671  {
672  AssertDimension(writable_rows.size(), row_parallel_partitioning.size());
674  reinit(row_parallel_partitioning.size(),
675  col_parallel_partitioning.size());
676  for (size_type i=0; i<row_parallel_partitioning.size(); ++i)
677  for (size_type j=0; j<col_parallel_partitioning.size(); ++j)
678  this->block(i,j).reinit(row_parallel_partitioning[i],
679  col_parallel_partitioning[j],
680  writable_rows[i],
681  communicator);
682  this->collect_sizes();
683  }
684 
685 }
686 
687 #endif
688 
691 #ifdef DEAL_II_WITH_TRILINOS
693 #endif
694 
695 DEAL_II_NAMESPACE_CLOSE
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
bool exists(const size_type i, const size_type j) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1052
void print(std::ostream &out) const
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:54
std::size_t memory_consumption() const
size_type size() const
Definition: index_set.h:1247
size_type n_nonzero_elements() const
SparsityPatternType & block(const size_type row, const size_type column)
void reinit(const size_type n_block_rows, const size_type n_block_columns)
types::global_dof_index size_type
void reinit(const size_type n_block_rows, const size_type n_block_columns)
const IndexSet & row_index_set() const
#define Assert(cond, exc)
Definition: exceptions.h:294
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
BlockSparsityPatternBase & operator=(const BlockSparsityPatternBase &)
size_type max_entries_per_row() const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
size_type block_size(const unsigned int i) const
void reinit(const size_type m, const size_type n, const unsigned int max_per_row)
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Table< 2, SmartPointer< SparsityPatternType, BlockSparsityPatternBase< SparsityPatternType > > > sub_objects
void print_gnuplot(std::ostream &out) const
bool is_element(const size_type index) const
Definition: index_set.h:1317
unsigned int size() const
void copy_from(const BlockDynamicSparsityPattern &dsp)
size_type local_to_global(const unsigned int block, const size_type index) const