Reference documentation for deal.II version Git e1b67c38dc 2020-08-13 12:19:15 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_sparsity_pattern.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2020 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.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
18 
19 #ifdef DEAL_II_WITH_TRILINOS
20 
21 # include <deal.II/base/mpi.h>
22 # include <deal.II/base/utilities.h>
23 
26 
27 # include <Epetra_Export.h>
28 
30 
31 namespace TrilinosWrappers
32 {
33  namespace SparsityPatternIterators
34  {
35  void
37  {
38  // if we are asked to visit the past-the-end line, then simply
39  // release all our caches and go on with life
40  if (this->a_row == sparsity_pattern->n_rows())
41  {
42  colnum_cache.reset();
43  return;
44  }
45 
46  // otherwise first flush Trilinos caches if necessary
49 
50  colnum_cache = std::make_shared<std::vector<size_type>>(
52 
53  if (colnum_cache->size() > 0)
54  {
55  // get a representation of the present row
56  int ncols;
57  const int ierr = sparsity_pattern->graph->ExtractGlobalRowCopy(
58  this->a_row,
59  colnum_cache->size(),
60  ncols,
61  reinterpret_cast<TrilinosWrappers::types::int_type *>(
62  const_cast<size_type *>(colnum_cache->data())));
63  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
64  AssertThrow(static_cast<std::vector<size_type>::size_type>(ncols) ==
65  colnum_cache->size(),
67  }
68  }
69  } // namespace SparsityPatternIterators
70 
71 
72  // The constructor is actually the
73  // only point where we have to check
74  // whether we build a serial or a
75  // parallel Trilinos matrix.
76  // Actually, it does not even matter
77  // how many threads there are, but
78  // only if we use an MPI compiler or
79  // a standard compiler. So, even one
80  // thread on a configuration with
81  // MPI will still get a parallel
82  // interface.
84  {
85  column_space_map =
86  std::make_unique<Epetra_Map>(TrilinosWrappers::types::int_type(0),
89  graph = std::make_unique<Epetra_FECrsGraph>(View,
90  *column_space_map,
91  *column_space_map,
92  0);
93  graph->FillComplete();
94  }
95 
96 
97 
99  const size_type n,
100  const size_type n_entries_per_row)
101  {
102  reinit(m, n, n_entries_per_row);
103  }
104 
105 
106 
108  const size_type m,
109  const size_type n,
110  const std::vector<size_type> &n_entries_per_row)
111  {
112  reinit(m, n, n_entries_per_row);
113  }
114 
115 
116 
118  : Subscriptor(std::move(other))
119  , column_space_map(std::move(other.column_space_map))
120  , graph(std::move(other.graph))
121  , nonlocal_graph(std::move(other.nonlocal_graph))
122  {}
123 
124 
125 
126  // Copy function only works if the
127  // sparsity pattern is empty.
129  : Subscriptor()
130  , column_space_map(new Epetra_Map(TrilinosWrappers::types::int_type(0),
132  Utilities::Trilinos::comm_self()))
133  , graph(
134  new Epetra_FECrsGraph(View, *column_space_map, *column_space_map, 0))
135  {
136  (void)input_sparsity;
137  Assert(input_sparsity.n_rows() == 0,
138  ExcMessage(
139  "Copy constructor only works for empty sparsity patterns."));
140  }
141 
142 
143 
144  SparsityPattern::SparsityPattern(const IndexSet &parallel_partitioning,
145  const MPI_Comm &communicator,
146  const size_type n_entries_per_row)
147  {
148  reinit(parallel_partitioning,
149  parallel_partitioning,
150  communicator,
151  n_entries_per_row);
152  }
153 
154 
155 
157  const IndexSet & parallel_partitioning,
158  const MPI_Comm & communicator,
159  const std::vector<size_type> &n_entries_per_row)
160  {
161  reinit(parallel_partitioning,
162  parallel_partitioning,
163  communicator,
164  n_entries_per_row);
165  }
166 
167 
168 
169  SparsityPattern::SparsityPattern(const IndexSet &row_parallel_partitioning,
170  const IndexSet &col_parallel_partitioning,
171  const MPI_Comm &communicator,
172  const size_type n_entries_per_row)
173  {
174  reinit(row_parallel_partitioning,
175  col_parallel_partitioning,
176  communicator,
177  n_entries_per_row);
178  }
179 
180 
181 
183  const IndexSet & row_parallel_partitioning,
184  const IndexSet & col_parallel_partitioning,
185  const MPI_Comm & communicator,
186  const std::vector<size_type> &n_entries_per_row)
187  {
188  reinit(row_parallel_partitioning,
189  col_parallel_partitioning,
190  communicator,
191  n_entries_per_row);
192  }
193 
194 
195 
196  SparsityPattern::SparsityPattern(const IndexSet &row_parallel_partitioning,
197  const IndexSet &col_parallel_partitioning,
198  const IndexSet &writable_rows,
199  const MPI_Comm &communicator,
200  const size_type n_max_entries_per_row)
201  {
202  reinit(row_parallel_partitioning,
203  col_parallel_partitioning,
204  writable_rows,
205  communicator,
206  n_max_entries_per_row);
207  }
208 
209 
210 
211  void
213  const size_type n,
214  const size_type n_entries_per_row)
215  {
218  MPI_COMM_SELF,
219  n_entries_per_row);
220  }
221 
222 
223 
224  void
226  const size_type n,
227  const std::vector<size_type> &n_entries_per_row)
228  {
231  MPI_COMM_SELF,
232  n_entries_per_row);
233  }
234 
235 
236 
237  namespace
238  {
240 
241  void
242  reinit_sp(const Epetra_Map & row_map,
243  const Epetra_Map & col_map,
244  const size_type n_entries_per_row,
245  std::unique_ptr<Epetra_Map> & column_space_map,
246  std::unique_ptr<Epetra_FECrsGraph> &graph,
247  std::unique_ptr<Epetra_CrsGraph> & nonlocal_graph)
248  {
249  Assert(row_map.IsOneToOne(),
250  ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
251  "the maps of different processors."));
252  Assert(col_map.IsOneToOne(),
253  ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
254  "the maps of different processors."));
255 
256  nonlocal_graph.reset();
257  graph.reset();
258  column_space_map = std::make_unique<Epetra_Map>(col_map);
259 
260  // for more than one processor, need to specify only row map first and
261  // let the matrix entries decide about the column map (which says which
262  // columns are present in the matrix, not to be confused with the
263  // col_map that tells how the domain dofs of the matrix will be
264  // distributed). for only one processor, we can directly assign the
265  // columns as well. If we use a recent Trilinos version, we can also
266  // require building a non-local graph which gives us thread-safe
267  // initialization.
268  if (row_map.Comm().NumProc() > 1)
269  graph = std::make_unique<Epetra_FECrsGraph>(
270  Copy, row_map, n_entries_per_row, false
271  // TODO: Check which new Trilinos version supports this...
272  // Remember to change tests/trilinos/assemble_matrix_parallel_07, too.
273  //#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
274  // , true
275  //#endif
276  );
277  else
278  graph = std::make_unique<Epetra_FECrsGraph>(
279  Copy, row_map, col_map, n_entries_per_row, false);
280  }
281 
282 
283 
284  void
285  reinit_sp(const Epetra_Map & row_map,
286  const Epetra_Map & col_map,
287  const std::vector<size_type> & n_entries_per_row,
288  std::unique_ptr<Epetra_Map> & column_space_map,
289  std::unique_ptr<Epetra_FECrsGraph> &graph,
290  std::unique_ptr<Epetra_CrsGraph> & nonlocal_graph)
291  {
292  Assert(row_map.IsOneToOne(),
293  ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
294  "the maps of different processors."));
295  Assert(col_map.IsOneToOne(),
296  ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
297  "the maps of different processors."));
298 
299  // release memory before reallocation
300  nonlocal_graph.reset();
301  graph.reset();
302  AssertDimension(n_entries_per_row.size(),
304 
305  column_space_map = std::make_unique<Epetra_Map>(col_map);
306  std::vector<int> local_entries_per_row(
309  for (unsigned int i = 0; i < local_entries_per_row.size(); ++i)
310  local_entries_per_row[i] =
311  n_entries_per_row[TrilinosWrappers::min_my_gid(row_map) + i];
312 
313  if (row_map.Comm().NumProc() > 1)
314  graph = std::make_unique<Epetra_FECrsGraph>(
315  Copy, row_map, local_entries_per_row.data(), false
316  // TODO: Check which new Trilinos version supports this...
317  // Remember to change tests/trilinos/assemble_matrix_parallel_07, too.
318  //#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
319  // , true
320  //#endif
321  );
322  else
323  graph = std::make_unique<Epetra_FECrsGraph>(
324  Copy, row_map, col_map, local_entries_per_row.data(), false);
325  }
326 
327 
328 
329  template <typename SparsityPatternType>
330  void
331  reinit_sp(const Epetra_Map & row_map,
332  const Epetra_Map & col_map,
333  const SparsityPatternType & sp,
334  const bool exchange_data,
335  std::unique_ptr<Epetra_Map> & column_space_map,
336  std::unique_ptr<Epetra_FECrsGraph> &graph,
337  std::unique_ptr<Epetra_CrsGraph> & nonlocal_graph)
338  {
339  nonlocal_graph.reset();
340  graph.reset();
341 
342  AssertDimension(sp.n_rows(),
344  AssertDimension(sp.n_cols(),
346 
347  column_space_map = std::make_unique<Epetra_Map>(col_map);
348 
349  Assert(row_map.LinearMap() == true,
350  ExcMessage(
351  "This function only works if the row map is contiguous."));
352 
353  const size_type first_row = TrilinosWrappers::min_my_gid(row_map),
354  last_row = TrilinosWrappers::max_my_gid(row_map) + 1;
355  std::vector<int> n_entries_per_row(last_row - first_row);
356 
357  // Trilinos wants the row length as an int this is hopefully never going
358  // to be a problem.
359  for (size_type row = first_row; row < last_row; ++row)
360  n_entries_per_row[row - first_row] =
361  static_cast<int>(sp.row_length(row));
362 
363  if (row_map.Comm().NumProc() > 1)
364  graph = std::make_unique<Epetra_FECrsGraph>(Copy,
365  row_map,
366  n_entries_per_row.data(),
367  false);
368  else
369  graph = std::make_unique<Epetra_FECrsGraph>(
370  Copy, row_map, col_map, n_entries_per_row.data(), false);
371 
372  AssertDimension(sp.n_rows(), n_global_rows(*graph));
373 
374  std::vector<TrilinosWrappers::types::int_type> row_indices;
375 
376  // Include possibility to exchange data since DynamicSparsityPattern is
377  // able to do so
378  if (exchange_data == false)
379  for (size_type row = first_row; row < last_row; ++row)
380  {
382  sp.row_length(row);
383  if (row_length == 0)
384  continue;
385 
386  row_indices.resize(row_length, -1);
387  {
388  typename SparsityPatternType::iterator p = sp.begin(row);
389  // avoid incrementing p over the end of the current row because
390  // it is slow for DynamicSparsityPattern in parallel
391  for (int col = 0; col < row_length;)
392  {
393  row_indices[col++] = p->column();
394  if (col < row_length)
395  ++p;
396  }
397  }
398  graph->Epetra_CrsGraph::InsertGlobalIndices(row,
399  row_length,
400  row_indices.data());
401  }
402  else
403  for (size_type row = 0; row < sp.n_rows(); ++row)
404  {
406  sp.row_length(row);
407  if (row_length == 0)
408  continue;
409 
410  row_indices.resize(row_length, -1);
411  {
412  typename SparsityPatternType::iterator p = sp.begin(row);
413  // avoid incrementing p over the end of the current row because
414  // it is slow for DynamicSparsityPattern in parallel
415  for (int col = 0; col < row_length;)
416  {
417  row_indices[col++] = p->column();
418  if (col < row_length)
419  ++p;
420  }
421  }
422  const TrilinosWrappers::types::int_type trilinos_row = row;
423  graph->InsertGlobalIndices(1,
424  &trilinos_row,
425  row_length,
426  row_indices.data());
427  }
428 
429  // TODO A dynamic_cast fails here, this is suspicious.
430  const auto &range_map =
431  static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
432  int ierr = graph->GlobalAssemble(*column_space_map, range_map, true);
433  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
434 
435  ierr = graph->OptimizeStorage();
436  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
437  }
438  } // namespace
439 
440 
441 
442  void
443  SparsityPattern::reinit(const IndexSet &parallel_partitioning,
444  const MPI_Comm &communicator,
445  const size_type n_entries_per_row)
446  {
447  Epetra_Map map =
448  parallel_partitioning.make_trilinos_map(communicator, false);
449  reinit_sp(
450  map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
451  }
452 
453 
454 
455  void
456  SparsityPattern::reinit(const IndexSet & parallel_partitioning,
457  const MPI_Comm & communicator,
458  const std::vector<size_type> &n_entries_per_row)
459  {
460  Epetra_Map map =
461  parallel_partitioning.make_trilinos_map(communicator, false);
462  reinit_sp(
463  map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
464  }
465 
466 
467 
468  void
469  SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
470  const IndexSet &col_parallel_partitioning,
471  const MPI_Comm &communicator,
472  const size_type n_entries_per_row)
473  {
474  Epetra_Map row_map =
475  row_parallel_partitioning.make_trilinos_map(communicator, false);
476  Epetra_Map col_map =
477  col_parallel_partitioning.make_trilinos_map(communicator, false);
478  reinit_sp(row_map,
479  col_map,
480  n_entries_per_row,
482  graph,
484  }
485 
486 
487 
488  void
489  SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
490  const IndexSet &col_parallel_partitioning,
491  const MPI_Comm &communicator,
492  const std::vector<size_type> &n_entries_per_row)
493  {
494  Epetra_Map row_map =
495  row_parallel_partitioning.make_trilinos_map(communicator, false);
496  Epetra_Map col_map =
497  col_parallel_partitioning.make_trilinos_map(communicator, false);
498  reinit_sp(row_map,
499  col_map,
500  n_entries_per_row,
502  graph,
504  }
505 
506 
507 
508  void
509  SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
510  const IndexSet &col_parallel_partitioning,
511  const IndexSet &writable_rows,
512  const MPI_Comm &communicator,
513  const size_type n_entries_per_row)
514  {
515  Epetra_Map row_map =
516  row_parallel_partitioning.make_trilinos_map(communicator, false);
517  Epetra_Map col_map =
518  col_parallel_partitioning.make_trilinos_map(communicator, false);
519  reinit_sp(row_map,
520  col_map,
521  n_entries_per_row,
523  graph,
525 
526  IndexSet nonlocal_partitioner = writable_rows;
527  AssertDimension(nonlocal_partitioner.size(),
528  row_parallel_partitioning.size());
529 # ifdef DEBUG
530  {
531  IndexSet tmp = writable_rows & row_parallel_partitioning;
532  Assert(tmp == row_parallel_partitioning,
533  ExcMessage(
534  "The set of writable rows passed to this method does not "
535  "contain the locally owned rows, which is not allowed."));
536  }
537 # endif
538  nonlocal_partitioner.subtract_set(row_parallel_partitioning);
539  if (Utilities::MPI::n_mpi_processes(communicator) > 1)
540  {
541  Epetra_Map nonlocal_map =
542  nonlocal_partitioner.make_trilinos_map(communicator, true);
544  std::make_unique<Epetra_CrsGraph>(Copy, nonlocal_map, 0);
545  }
546  else
547  Assert(nonlocal_partitioner.n_elements() == 0, ExcInternalError());
548  }
549 
550 
551 
552  template <typename SparsityPatternType>
553  void
555  const IndexSet & row_parallel_partitioning,
556  const IndexSet & col_parallel_partitioning,
557  const SparsityPatternType &nontrilinos_sparsity_pattern,
558  const MPI_Comm & communicator,
559  const bool exchange_data)
560  {
561  Epetra_Map row_map =
562  row_parallel_partitioning.make_trilinos_map(communicator, false);
563  Epetra_Map col_map =
564  col_parallel_partitioning.make_trilinos_map(communicator, false);
565  reinit_sp(row_map,
566  col_map,
567  nontrilinos_sparsity_pattern,
568  exchange_data,
570  graph,
572  }
573 
574 
575 
576  template <typename SparsityPatternType>
577  void
579  const IndexSet & parallel_partitioning,
580  const SparsityPatternType &nontrilinos_sparsity_pattern,
581  const MPI_Comm & communicator,
582  const bool exchange_data)
583  {
584  Epetra_Map map =
585  parallel_partitioning.make_trilinos_map(communicator, false);
586  reinit_sp(map,
587  map,
588  nontrilinos_sparsity_pattern,
589  exchange_data,
591  graph,
593  }
594 
595 
596 
599  {
600  Assert(false, ExcNotImplemented());
601  return *this;
602  }
603 
604 
605 
606  void
608  {
609  column_space_map = std::make_unique<Epetra_Map>(*sp.column_space_map);
610  graph = std::make_unique<Epetra_FECrsGraph>(*sp.graph);
611 
612  if (sp.nonlocal_graph.get() != nullptr)
613  nonlocal_graph = std::make_unique<Epetra_CrsGraph>(*sp.nonlocal_graph);
614  else
615  nonlocal_graph.reset();
616  }
617 
618 
619 
620  template <typename SparsityPatternType>
621  void
622  SparsityPattern::copy_from(const SparsityPatternType &sp)
623  {
624  const Epetra_Map rows(TrilinosWrappers::types::int_type(sp.n_rows()),
625  0,
627  const Epetra_Map columns(TrilinosWrappers::types::int_type(sp.n_cols()),
628  0,
630 
631  reinit_sp(
632  rows, columns, sp, false, column_space_map, graph, nonlocal_graph);
633  }
634 
635 
636 
637  void
639  {
640  // When we clear the matrix, reset
641  // the pointer and generate an
642  // empty sparsity pattern.
644  std::make_unique<Epetra_Map>(TrilinosWrappers::types::int_type(0),
647  graph = std::make_unique<Epetra_FECrsGraph>(View,
650  0);
651  graph->FillComplete();
652 
653  nonlocal_graph.reset();
654  }
655 
656 
657 
658  void
660  {
661  int ierr;
663  if (nonlocal_graph.get() != nullptr)
664  {
665  if (nonlocal_graph->IndicesAreGlobal() == false &&
666  nonlocal_graph->RowMap().NumMyElements() > 0 &&
668  {
669  // Insert dummy element at (row, column) that corresponds to row 0
670  // in local index counting.
674 
675  // in case we have a square sparsity pattern, add the entry on the
676  // diagonal
679  column = row;
680  // if not, take a column index that we have ourselves since we
681  // know for sure it is there (and it will not create spurious
682  // messages to many ranks like putting index 0 on many processors)
683  else if (column_space_map->NumMyElements() > 0)
685  ierr = nonlocal_graph->InsertGlobalIndices(row, 1, &column);
686  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
687  }
688  Assert(nonlocal_graph->RowMap().NumMyElements() == 0 ||
690  nonlocal_graph->IndicesAreGlobal() == true,
691  ExcInternalError());
692 
693  ierr =
694  nonlocal_graph->FillComplete(*column_space_map, graph->RangeMap());
695  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
696  ierr = nonlocal_graph->OptimizeStorage();
697  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
698  Epetra_Export exporter(nonlocal_graph->RowMap(), graph->RowMap());
699  ierr = graph->Export(*nonlocal_graph, exporter, Add);
700  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
701  ierr = graph->FillComplete(*column_space_map, graph->RangeMap());
702  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
703  }
704  else
705  {
706  // TODO A dynamic_cast fails here, this is suspicious.
707  const auto &range_map =
708  static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
709  ierr = graph->GlobalAssemble(*column_space_map, range_map, true);
710  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
711  }
712 
713  ierr = graph->OptimizeStorage();
714  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
715  }
716 
717 
718 
719  bool
721  {
722  return graph->RowMap().LID(
723  static_cast<TrilinosWrappers::types::int_type>(i)) != -1;
724  }
725 
726 
727 
728  bool
730  {
731  if (!row_is_stored_locally(i))
732  {
733  return false;
734  }
735  else
736  {
737  // Extract local indices in
738  // the matrix.
739  int trilinos_i =
740  graph->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
741  trilinos_j =
742  graph->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
743 
744  // Check whether the matrix
745  // already is transformed to
746  // local indices.
747  if (graph->Filled() == false)
748  {
749  int nnz_present = graph->NumGlobalIndices(i);
750  int nnz_extracted;
752 
753  // Generate the view and make
754  // sure that we have not generated
755  // an error.
756  // TODO: trilinos_i is the local row index -> it is an int but
757  // ExtractGlobalRowView requires trilinos_i to be the global row
758  // index and thus it should be a long long int
759  int ierr = graph->ExtractGlobalRowView(trilinos_i,
760  nnz_extracted,
761  col_indices);
762  (void)ierr;
763  Assert(ierr == 0, ExcTrilinosError(ierr));
764  Assert(nnz_present == nnz_extracted,
765  ExcDimensionMismatch(nnz_present, nnz_extracted));
766 
767  // Search the index
768  const std::ptrdiff_t local_col_index =
769  std::find(col_indices, col_indices + nnz_present, trilinos_j) -
770  col_indices;
771 
772  if (local_col_index == nnz_present)
773  return false;
774  }
775  else
776  {
777  // Prepare pointers for extraction
778  // of a view of the row.
779  int nnz_present = graph->NumGlobalIndices(i);
780  int nnz_extracted;
781  int *col_indices;
782 
783  // Generate the view and make
784  // sure that we have not generated
785  // an error.
786  int ierr =
787  graph->ExtractMyRowView(trilinos_i, nnz_extracted, col_indices);
788  (void)ierr;
789  Assert(ierr == 0, ExcTrilinosError(ierr));
790 
791  Assert(nnz_present == nnz_extracted,
792  ExcDimensionMismatch(nnz_present, nnz_extracted));
793 
794  // Search the index
795  const std::ptrdiff_t local_col_index =
796  std::find(col_indices, col_indices + nnz_present, trilinos_j) -
797  col_indices;
798 
799  if (local_col_index == nnz_present)
800  return false;
801  }
802  }
803 
804  return true;
805  }
806 
807 
808 
811  {
812  size_type local_b = 0;
814  for (int i = 0; i < static_cast<int>(local_size()); ++i)
815  {
816  int *indices;
817  int num_entries;
818  graph->ExtractMyRowView(i, num_entries, indices);
819  for (unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
820  ++j)
821  {
822  if (static_cast<size_type>(std::abs(i - indices[j])) > local_b)
823  local_b = std::abs(i - indices[j]);
824  }
825  }
826  graph->Comm().MaxAll(reinterpret_cast<TrilinosWrappers::types::int_type *>(
827  &local_b),
828  &global_b,
829  1);
830  return static_cast<size_type>(global_b);
831  }
832 
833 
834 
837  {
839  return n_rows;
840  }
841 
842 
843 
846  {
848  if (graph->Filled() == true)
849  n_cols = n_global_cols(*graph);
850  else
852 
853  return n_cols;
854  }
855 
856 
857 
858  unsigned int
860  {
861  int n_rows = graph->NumMyRows();
862 
863  return n_rows;
864  }
865 
866 
867 
868  std::pair<SparsityPattern::size_type, SparsityPattern::size_type>
870  {
872  begin = TrilinosWrappers::min_my_gid(graph->RowMap());
873  end = TrilinosWrappers::max_my_gid(graph->RowMap()) + 1;
874 
875  return std::make_pair(begin, end);
876  }
877 
878 
879 
882  {
884 
885  return static_cast<size_type>(nnz);
886  }
887 
888 
889 
890  unsigned int
892  {
893  int nnz = graph->MaxNumIndices();
894 
895  return static_cast<unsigned int>(nnz);
896  }
897 
898 
899 
902  {
903  Assert(row < n_rows(), ExcInternalError());
904 
905  // Get a representation of the where the present row is located on
906  // the current processor
908  graph->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
909 
910  // On the processor who owns this row, we'll have a non-negative
911  // value for `local_row` and can ask for the length of the row.
912  if (local_row >= 0)
913  return graph->NumMyIndices(local_row);
914  else
915  return static_cast<size_type>(-1);
916  }
917 
918 
919 
920  const Epetra_Map &
922  {
923  // TODO A dynamic_cast fails here, this is suspicious.
924  const auto &domain_map =
925  static_cast<const Epetra_Map &>(graph->DomainMap()); // NOLINT
926  return domain_map;
927  }
928 
929 
930 
931  const Epetra_Map &
933  {
934  // TODO A dynamic_cast fails here, this is suspicious.
935  const auto &range_map =
936  static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
937  return range_map;
938  }
939 
940 
941 
942  MPI_Comm
944  {
945 # ifdef DEAL_II_WITH_MPI
946 
947  const Epetra_MpiComm *mpi_comm =
948  dynamic_cast<const Epetra_MpiComm *>(&graph->RangeMap().Comm());
949  Assert(mpi_comm != nullptr, ExcInternalError());
950  return mpi_comm->Comm();
951 # else
952 
953  return MPI_COMM_SELF;
954 
955 # endif
956  }
957 
958 
959 
960  void
962  {
963  Assert(false, ExcNotImplemented());
964  }
965 
966 
967 
968  // As of now, no particularly neat
969  // output is generated in case of
970  // multiple processors.
971  void
972  SparsityPattern::print(std::ostream &out,
973  const bool write_extended_trilinos_info) const
974  {
975  if (write_extended_trilinos_info)
976  out << *graph;
977  else
978  {
979  int *indices;
980  int num_entries;
981 
982  for (int i = 0; i < graph->NumMyRows(); ++i)
983  {
984  graph->ExtractMyRowView(i, num_entries, indices);
985  for (int j = 0; j < num_entries; ++j)
986  out << "(" << TrilinosWrappers::global_index(graph->RowMap(), i)
987  << ","
988  << TrilinosWrappers::global_index(graph->ColMap(), indices[j])
989  << ") " << std::endl;
990  }
991  }
992 
993  AssertThrow(out, ExcIO());
994  }
995 
996 
997 
998  void
999  SparsityPattern::print_gnuplot(std::ostream &out) const
1000  {
1001  Assert(graph->Filled() == true, ExcInternalError());
1002  for (::types::global_dof_index row = 0; row < local_size(); ++row)
1003  {
1004  int *indices;
1005  int num_entries;
1006  graph->ExtractMyRowView(row, num_entries, indices);
1007 
1008  Assert(num_entries >= 0, ExcInternalError());
1009  // avoid sign comparison warning
1010  const ::types::global_dof_index num_entries_ = num_entries;
1011  for (::types::global_dof_index j = 0; j < num_entries_; ++j)
1012  // while matrix entries are usually
1013  // written (i,j), with i vertical and
1014  // j horizontal, gnuplot output is
1015  // x-y, that is we have to exchange
1016  // the order of output
1017  out << static_cast<int>(
1018  TrilinosWrappers::global_index(graph->ColMap(), indices[j]))
1019  << " "
1020  << -static_cast<int>(
1021  TrilinosWrappers::global_index(graph->RowMap(), row))
1022  << std::endl;
1023  }
1024 
1025  AssertThrow(out, ExcIO());
1026  }
1027 
1028  // TODO: Implement!
1029  std::size_t
1031  {
1032  Assert(false, ExcNotImplemented());
1033  return 0;
1034  }
1035 
1036 
1037 # ifndef DOXYGEN
1038  // explicit instantiations
1039  //
1040  template void
1041  SparsityPattern::copy_from(const ::SparsityPattern &);
1042  template void
1043  SparsityPattern::copy_from(const ::DynamicSparsityPattern &);
1044 
1045  template void
1047  const ::SparsityPattern &,
1048  const MPI_Comm &,
1049  bool);
1050  template void
1052  const ::DynamicSparsityPattern &,
1053  const MPI_Comm &,
1054  bool);
1055 
1056 
1057  template void
1059  const IndexSet &,
1060  const ::SparsityPattern &,
1061  const MPI_Comm &,
1062  bool);
1063  template void
1065  const IndexSet &,
1066  const ::DynamicSparsityPattern &,
1067  const MPI_Comm &,
1068  bool);
1069 # endif
1070 
1071 } // namespace TrilinosWrappers
1072 
1074 
1075 #endif // DEAL_II_WITH_TRILINOS
static ::ExceptionBase & ExcTrilinosError(int arg1)
const Epetra_Map & domain_partitioner() const
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
typename ::internal::FEValuesViews::ViewType< dim, spacedim, Extractor >::type View
Definition: fe_values.h:1970
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
types::global_dof_index size_type
Definition: cuda_kernels.h:45
std::shared_ptr< const std::vector< size_type > > colnum_cache
static ::ExceptionBase & ExcIO()
size_type row_length(const size_type row) const
TrilinosWrappers::types::int_type min_my_gid(const Epetra_BlockMap &map)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1521
const Epetra_Comm & comm_self()
Definition: utilities.cc:1114
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
const_iterator begin() const
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:596
size_type size() const
Definition: index_set.h:1632
bool exists(const size_type i, const size_type j) const
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::types::int_type n_global_rows(const Epetra_CrsGraph &graph)
Definition: types.h:31
std::unique_ptr< Epetra_Map > column_space_map
void subtract_set(const IndexSet &other)
Definition: index_set.cc:258
#define Assert(cond, exc)
Definition: exceptions.h:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
TrilinosWrappers::types::int_type n_global_cols(const Epetra_CrsGraph &graph)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
void print_gnuplot(std::ostream &out) const
const_iterator end() const
TrilinosWrappers::types::int_type n_global_entries(const Epetra_CrsGraph &graph)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
TrilinosWrappers::types::int_type n_global_elements(const Epetra_BlockMap &map)
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1011
unsigned int global_dof_index
Definition: types.h:76
SparsityPattern & operator=(const SparsityPattern &input_sparsity_pattern)
Definition: cuda.h:31
void copy_from(const SparsityPattern &input_sparsity_pattern)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
TrilinosWrappers::types::int_type max_my_gid(const Epetra_BlockMap &map)
static ::ExceptionBase & ExcNotImplemented()
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
void reinit(const size_type m, const size_type n, const size_type n_entries_per_row=0)
const Epetra_Map & range_partitioner() const
size_type n_elements() const
Definition: index_set.h:1830
std::pair< size_type, size_type > local_range() const
std::unique_ptr< Epetra_FECrsGraph > graph
static ::ExceptionBase & ExcInternalError()
bool row_is_stored_locally(const size_type i) const