Reference documentation for deal.II version GIT e22cb6a53e 2023-09-21 14:00:02+00:00
\(\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 - 2023 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>
23 
26 
27 # include <Epetra_Export.h>
28 
29 # include <limits>
30 
32 
33 namespace TrilinosWrappers
34 {
35  namespace SparsityPatternIterators
36  {
37  void
39  {
40  // if we are asked to visit the past-the-end line, then simply
41  // release all our caches and go on with life
42  if (this->a_row == sparsity_pattern->n_rows())
43  {
44  colnum_cache.reset();
45  return;
46  }
47 
48  // otherwise first flush Trilinos caches if necessary
51 
52  colnum_cache = std::make_shared<std::vector<size_type>>(
54 
55  if (colnum_cache->size() > 0)
56  {
57  // get a representation of the present row
58  int ncols;
59  const int ierr = sparsity_pattern->graph->ExtractGlobalRowCopy(
60  this->a_row,
61  colnum_cache->size(),
62  ncols,
63  reinterpret_cast<TrilinosWrappers::types::int_type *>(
64  const_cast<size_type *>(colnum_cache->data())));
65  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
66  AssertThrow(static_cast<std::vector<size_type>::size_type>(ncols) ==
67  colnum_cache->size(),
69  }
70  }
71  } // namespace SparsityPatternIterators
72 
73 
74  // The constructor is actually the
75  // only point where we have to check
76  // whether we build a serial or a
77  // parallel Trilinos matrix.
78  // Actually, it does not even matter
79  // how many threads there are, but
80  // only if we use an MPI compiler or
81  // a standard compiler. So, even one
82  // thread on a configuration with
83  // MPI will still get a parallel
84  // interface.
86  {
88  std::make_unique<Epetra_Map>(TrilinosWrappers::types::int_type(0),
91  graph = std::make_unique<Epetra_FECrsGraph>(View,
94  0);
95  graph->FillComplete();
96  }
97 
98 
99 
101  const size_type n,
102  const size_type n_entries_per_row)
103  {
104  reinit(m, n, n_entries_per_row);
105  }
106 
107 
108 
110  const size_type m,
111  const size_type n,
112  const std::vector<size_type> &n_entries_per_row)
113  {
114  reinit(m, n, n_entries_per_row);
115  }
116 
117 
118 
120  : SparsityPatternBase(std::move(other))
121  , column_space_map(std::move(other.column_space_map))
122  , graph(std::move(other.graph))
123  , nonlocal_graph(std::move(other.nonlocal_graph))
124  {}
125 
126 
127 
128  // Copy function only works if the
129  // sparsity pattern is empty.
131  : SparsityPatternBase(input_sparsity)
132  , column_space_map(new Epetra_Map(TrilinosWrappers::types::int_type(0),
134  Utilities::Trilinos::comm_self()))
135  , graph(
136  new Epetra_FECrsGraph(View, *column_space_map, *column_space_map, 0))
137  {
138  (void)input_sparsity;
139  Assert(input_sparsity.n_rows() == 0,
140  ExcMessage(
141  "Copy constructor only works for empty sparsity patterns."));
142  }
143 
144 
145 
146  SparsityPattern::SparsityPattern(const IndexSet &parallel_partitioning,
147  const MPI_Comm communicator,
148  const size_type n_entries_per_row)
149  {
150  reinit(parallel_partitioning,
151  parallel_partitioning,
152  communicator,
153  n_entries_per_row);
154  }
155 
156 
157 
159  const IndexSet &parallel_partitioning,
160  const MPI_Comm communicator,
161  const std::vector<size_type> &n_entries_per_row)
162  {
163  reinit(parallel_partitioning,
164  parallel_partitioning,
165  communicator,
166  n_entries_per_row);
167  }
168 
169 
170 
171  SparsityPattern::SparsityPattern(const IndexSet &row_parallel_partitioning,
172  const IndexSet &col_parallel_partitioning,
173  const MPI_Comm communicator,
174  const size_type n_entries_per_row)
175  {
176  reinit(row_parallel_partitioning,
177  col_parallel_partitioning,
178  communicator,
179  n_entries_per_row);
180  }
181 
182 
183 
185  const IndexSet &row_parallel_partitioning,
186  const IndexSet &col_parallel_partitioning,
187  const MPI_Comm communicator,
188  const std::vector<size_type> &n_entries_per_row)
189  {
190  reinit(row_parallel_partitioning,
191  col_parallel_partitioning,
192  communicator,
193  n_entries_per_row);
194  }
195 
196 
197 
198  SparsityPattern::SparsityPattern(const IndexSet &row_parallel_partitioning,
199  const IndexSet &col_parallel_partitioning,
200  const IndexSet &writable_rows,
201  const MPI_Comm communicator,
202  const size_type n_max_entries_per_row)
203  {
204  reinit(row_parallel_partitioning,
205  col_parallel_partitioning,
206  writable_rows,
207  communicator,
208  n_max_entries_per_row);
209  }
210 
211 
212 
213  void
215  const size_type n,
216  const size_type n_entries_per_row)
217  {
220  MPI_COMM_SELF,
221  n_entries_per_row);
222  }
223 
224 
225 
226  void
228  const size_type n,
229  const std::vector<size_type> &n_entries_per_row)
230  {
233  MPI_COMM_SELF,
234  n_entries_per_row);
235  }
236 
237 
238 
239  namespace
240  {
242 
243  void
244  reinit_sp(const Epetra_Map &row_map,
245  const Epetra_Map &col_map,
246  const size_type n_entries_per_row,
247  std::unique_ptr<Epetra_Map> &column_space_map,
248  std::unique_ptr<Epetra_FECrsGraph> &graph,
249  std::unique_ptr<Epetra_CrsGraph> &nonlocal_graph)
250  {
251  Assert(row_map.IsOneToOne(),
252  ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
253  "the maps of different processors."));
254  Assert(col_map.IsOneToOne(),
255  ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
256  "the maps of different processors."));
257 
258  nonlocal_graph.reset();
259  graph.reset();
260  column_space_map = std::make_unique<Epetra_Map>(col_map);
261 
263  TrilinosWrappers::min_my_gid(row_map)) *
264  std::uint64_t(n_entries_per_row) <
265  static_cast<std::uint64_t>(std::numeric_limits<int>::max()),
266  ExcMessage("The TrilinosWrappers use Epetra internally which "
267  "uses 'signed int' to represent local indices. "
268  "Therefore, only 2,147,483,647 nonzero matrix "
269  "entries can be stored on a single process, "
270  "but you are requesting more than that. "
271  "If possible, use more MPI processes."));
272 
273 
274  // for more than one processor, need to specify only row map first and
275  // let the matrix entries decide about the column map (which says which
276  // columns are present in the matrix, not to be confused with the
277  // col_map that tells how the domain dofs of the matrix will be
278  // distributed). for only one processor, we can directly assign the
279  // columns as well. If we use a recent Trilinos version, we can also
280  // require building a non-local graph which gives us thread-safe
281  // initialization.
282  if (row_map.Comm().NumProc() > 1)
283  graph = std::make_unique<Epetra_FECrsGraph>(
284  Copy, row_map, n_entries_per_row, false
285  // TODO: Check which new Trilinos version supports this...
286  // Remember to change tests/trilinos/assemble_matrix_parallel_07, too.
287  // #if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
288  // , true
289  // #endif
290  );
291  else
292  graph = std::make_unique<Epetra_FECrsGraph>(
293  Copy, row_map, col_map, n_entries_per_row, false);
294  }
295 
296 
297 
298  void
299  reinit_sp(const Epetra_Map &row_map,
300  const Epetra_Map &col_map,
301  const std::vector<size_type> &n_entries_per_row,
302  std::unique_ptr<Epetra_Map> &column_space_map,
303  std::unique_ptr<Epetra_FECrsGraph> &graph,
304  std::unique_ptr<Epetra_CrsGraph> &nonlocal_graph)
305  {
306  Assert(row_map.IsOneToOne(),
307  ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
308  "the maps of different processors."));
309  Assert(col_map.IsOneToOne(),
310  ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
311  "the maps of different processors."));
312 
313  // release memory before reallocation
314  nonlocal_graph.reset();
315  graph.reset();
316  AssertDimension(n_entries_per_row.size(),
318 
319  column_space_map = std::make_unique<Epetra_Map>(col_map);
320 
321  // Translate the vector of row lengths into one that only stores
322  // those entries that related to the locally stored rows of the matrix:
323  std::vector<int> local_entries_per_row(
326  for (unsigned int i = 0; i < local_entries_per_row.size(); ++i)
327  local_entries_per_row[i] =
328  n_entries_per_row[TrilinosWrappers::min_my_gid(row_map) + i];
329 
330  AssertThrow(std::accumulate(local_entries_per_row.begin(),
331  local_entries_per_row.end(),
332  std::uint64_t(0)) <
333  static_cast<std::uint64_t>(std::numeric_limits<int>::max()),
334  ExcMessage("The TrilinosWrappers use Epetra internally which "
335  "uses 'signed int' to represent local indices. "
336  "Therefore, only 2,147,483,647 nonzero matrix "
337  "entries can be stored on a single process, "
338  "but you are requesting more than that. "
339  "If possible, use more MPI processes."));
340 
341  if (row_map.Comm().NumProc() > 1)
342  graph = std::make_unique<Epetra_FECrsGraph>(
343  Copy, row_map, local_entries_per_row.data(), false
344  // TODO: Check which new Trilinos version supports this...
345  // Remember to change tests/trilinos/assemble_matrix_parallel_07, too.
346  // #if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
347  // , true
348  // #endif
349  );
350  else
351  graph = std::make_unique<Epetra_FECrsGraph>(
352  Copy, row_map, col_map, local_entries_per_row.data(), false);
353  }
354 
355 
356 
357  template <typename SparsityPatternType>
358  void
359  reinit_sp(const Epetra_Map &row_map,
360  const Epetra_Map &col_map,
361  const SparsityPatternType &sp,
362  const bool exchange_data,
363  std::unique_ptr<Epetra_Map> &column_space_map,
364  std::unique_ptr<Epetra_FECrsGraph> &graph,
365  std::unique_ptr<Epetra_CrsGraph> &nonlocal_graph)
366  {
367  nonlocal_graph.reset();
368  graph.reset();
369 
370  AssertDimension(sp.n_rows(),
372  AssertDimension(sp.n_cols(),
374 
375  column_space_map = std::make_unique<Epetra_Map>(col_map);
376 
377  Assert(row_map.LinearMap() == true,
378  ExcMessage(
379  "This function only works if the row map is contiguous."));
380 
381  const size_type first_row = TrilinosWrappers::min_my_gid(row_map),
382  last_row = TrilinosWrappers::max_my_gid(row_map) + 1;
383  std::vector<int> n_entries_per_row(last_row - first_row);
384 
385  // Trilinos wants the row length as an int this is hopefully never going
386  // to be a problem.
387  for (size_type row = first_row; row < last_row; ++row)
388  n_entries_per_row[row - first_row] =
389  static_cast<int>(sp.row_length(row));
390 
391  AssertThrow(std::accumulate(n_entries_per_row.begin(),
392  n_entries_per_row.end(),
393  std::uint64_t(0)) <
394  static_cast<std::uint64_t>(std::numeric_limits<int>::max()),
395  ExcMessage("The TrilinosWrappers use Epetra internally which "
396  "uses 'signed int' to represent local indices. "
397  "Therefore, only 2,147,483,647 nonzero matrix "
398  "entries can be stored on a single process, "
399  "but you are requesting more than that. "
400  "If possible, use more MPI processes."));
401 
402  if (row_map.Comm().NumProc() > 1)
403  graph = std::make_unique<Epetra_FECrsGraph>(Copy,
404  row_map,
405  n_entries_per_row.data(),
406  false);
407  else
408  graph = std::make_unique<Epetra_FECrsGraph>(
409  Copy, row_map, col_map, n_entries_per_row.data(), false);
410 
411  AssertDimension(sp.n_rows(), n_global_rows(*graph));
412 
413  std::vector<TrilinosWrappers::types::int_type> row_indices;
414 
415  // Include possibility to exchange data since DynamicSparsityPattern is
416  // able to do so
417  if (exchange_data == false)
418  for (size_type row = first_row; row < last_row; ++row)
419  {
420  const TrilinosWrappers::types::int_type row_length =
421  sp.row_length(row);
422  if (row_length == 0)
423  continue;
424 
425  row_indices.resize(row_length, -1);
426  {
427  typename SparsityPatternType::iterator p = sp.begin(row);
428  // avoid incrementing p over the end of the current row because
429  // it is slow for DynamicSparsityPattern in parallel
430  for (int col = 0; col < row_length;)
431  {
432  row_indices[col++] = p->column();
433  if (col < row_length)
434  ++p;
435  }
436  }
437  graph->Epetra_CrsGraph::InsertGlobalIndices(row,
438  row_length,
439  row_indices.data());
440  }
441  else
442  for (size_type row = 0; row < sp.n_rows(); ++row)
443  {
444  const TrilinosWrappers::types::int_type row_length =
445  sp.row_length(row);
446  if (row_length == 0)
447  continue;
448 
449  row_indices.resize(row_length, -1);
450  {
451  typename SparsityPatternType::iterator p = sp.begin(row);
452  // avoid incrementing p over the end of the current row because
453  // it is slow for DynamicSparsityPattern in parallel
454  for (int col = 0; col < row_length;)
455  {
456  row_indices[col++] = p->column();
457  if (col < row_length)
458  ++p;
459  }
460  }
461  const TrilinosWrappers::types::int_type trilinos_row = row;
462  graph->InsertGlobalIndices(1,
463  &trilinos_row,
464  row_length,
465  row_indices.data());
466  }
467 
468  // TODO A dynamic_cast fails here, this is suspicious.
469  const auto &range_map =
470  static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
471  int ierr = graph->GlobalAssemble(*column_space_map, range_map, true);
472  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
473 
474  ierr = graph->OptimizeStorage();
475  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
476  }
477  } // namespace
478 
479 
480 
481  void
482  SparsityPattern::reinit(const IndexSet &parallel_partitioning,
483  const MPI_Comm communicator,
484  const size_type n_entries_per_row)
485  {
486  SparsityPatternBase::resize(parallel_partitioning.size(),
487  parallel_partitioning.size());
488  Epetra_Map map =
489  parallel_partitioning.make_trilinos_map(communicator, false);
490  reinit_sp(
491  map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
492  }
493 
494 
495 
496  void
497  SparsityPattern::reinit(const IndexSet &parallel_partitioning,
498  const MPI_Comm communicator,
499  const std::vector<size_type> &n_entries_per_row)
500  {
501  SparsityPatternBase::resize(parallel_partitioning.size(),
502  parallel_partitioning.size());
503  Epetra_Map map =
504  parallel_partitioning.make_trilinos_map(communicator, false);
505  reinit_sp(
506  map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
507  }
508 
509 
510 
511  void
512  SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
513  const IndexSet &col_parallel_partitioning,
514  const MPI_Comm communicator,
515  const size_type n_entries_per_row)
516  {
517  SparsityPatternBase::resize(row_parallel_partitioning.size(),
518  col_parallel_partitioning.size());
519  Epetra_Map row_map =
520  row_parallel_partitioning.make_trilinos_map(communicator, false);
521  Epetra_Map col_map =
522  col_parallel_partitioning.make_trilinos_map(communicator, false);
523  reinit_sp(row_map,
524  col_map,
525  n_entries_per_row,
527  graph,
529  }
530 
531 
532 
533  void
534  SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
535  const IndexSet &col_parallel_partitioning,
536  const MPI_Comm communicator,
537  const std::vector<size_type> &n_entries_per_row)
538  {
539  SparsityPatternBase::resize(row_parallel_partitioning.size(),
540  col_parallel_partitioning.size());
541  Epetra_Map row_map =
542  row_parallel_partitioning.make_trilinos_map(communicator, false);
543  Epetra_Map col_map =
544  col_parallel_partitioning.make_trilinos_map(communicator, false);
545  reinit_sp(row_map,
546  col_map,
547  n_entries_per_row,
549  graph,
551  }
552 
553 
554 
555  void
556  SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
557  const IndexSet &col_parallel_partitioning,
558  const IndexSet &writable_rows,
559  const MPI_Comm communicator,
560  const size_type n_entries_per_row)
561  {
562  SparsityPatternBase::resize(row_parallel_partitioning.size(),
563  col_parallel_partitioning.size());
564  Epetra_Map row_map =
565  row_parallel_partitioning.make_trilinos_map(communicator, false);
566  Epetra_Map col_map =
567  col_parallel_partitioning.make_trilinos_map(communicator, false);
568  reinit_sp(row_map,
569  col_map,
570  n_entries_per_row,
572  graph,
574 
575  IndexSet nonlocal_partitioner = writable_rows;
576  AssertDimension(nonlocal_partitioner.size(),
577  row_parallel_partitioning.size());
578 # ifdef DEBUG
579  {
580  IndexSet tmp = writable_rows & row_parallel_partitioning;
581  Assert(tmp == row_parallel_partitioning,
582  ExcMessage(
583  "The set of writable rows passed to this method does not "
584  "contain the locally owned rows, which is not allowed."));
585  }
586 # endif
587  nonlocal_partitioner.subtract_set(row_parallel_partitioning);
588  if (Utilities::MPI::n_mpi_processes(communicator) > 1)
589  {
590  Epetra_Map nonlocal_map =
591  nonlocal_partitioner.make_trilinos_map(communicator, true);
593  std::make_unique<Epetra_CrsGraph>(Copy, nonlocal_map, 0);
594  }
595  else
596  Assert(nonlocal_partitioner.n_elements() == 0, ExcInternalError());
597  }
598 
599 
600 
601  template <typename SparsityPatternType>
602  void
604  const IndexSet &row_parallel_partitioning,
605  const IndexSet &col_parallel_partitioning,
606  const SparsityPatternType &nontrilinos_sparsity_pattern,
607  const MPI_Comm communicator,
608  const bool exchange_data)
609  {
610  SparsityPatternBase::resize(row_parallel_partitioning.size(),
611  col_parallel_partitioning.size());
612  Epetra_Map row_map =
613  row_parallel_partitioning.make_trilinos_map(communicator, false);
614  Epetra_Map col_map =
615  col_parallel_partitioning.make_trilinos_map(communicator, false);
616  reinit_sp(row_map,
617  col_map,
618  nontrilinos_sparsity_pattern,
619  exchange_data,
621  graph,
623  }
624 
625 
626 
627  template <typename SparsityPatternType>
628  void
630  const IndexSet &parallel_partitioning,
631  const SparsityPatternType &nontrilinos_sparsity_pattern,
632  const MPI_Comm communicator,
633  const bool exchange_data)
634  {
635  AssertDimension(nontrilinos_sparsity_pattern.n_rows(),
636  parallel_partitioning.size());
637  AssertDimension(nontrilinos_sparsity_pattern.n_cols(),
638  parallel_partitioning.size());
639  SparsityPatternBase::resize(parallel_partitioning.size(),
640  parallel_partitioning.size());
641  Epetra_Map map =
642  parallel_partitioning.make_trilinos_map(communicator, false);
643  reinit_sp(map,
644  map,
645  nontrilinos_sparsity_pattern,
646  exchange_data,
648  graph,
650  }
651 
652 
653 
656  {
657  Assert(false, ExcNotImplemented());
658  return *this;
659  }
660 
661 
662 
663  void
665  {
667  column_space_map = std::make_unique<Epetra_Map>(*sp.column_space_map);
668  graph = std::make_unique<Epetra_FECrsGraph>(*sp.graph);
669 
670  if (sp.nonlocal_graph.get() != nullptr)
671  nonlocal_graph = std::make_unique<Epetra_CrsGraph>(*sp.nonlocal_graph);
672  else
673  nonlocal_graph.reset();
674  }
675 
676 
677 
678  template <typename SparsityPatternType>
679  void
680  SparsityPattern::copy_from(const SparsityPatternType &sp)
681  {
682  SparsityPatternBase::resize(sp.n_rows(), sp.n_cols());
683  const Epetra_Map rows(TrilinosWrappers::types::int_type(sp.n_rows()),
684  0,
686  const Epetra_Map columns(TrilinosWrappers::types::int_type(sp.n_cols()),
687  0,
689 
690  reinit_sp(
691  rows, columns, sp, false, column_space_map, graph, nonlocal_graph);
692  }
693 
694 
695 
696  void
698  {
700  // When we clear the matrix, reset
701  // the pointer and generate an
702  // empty sparsity pattern.
704  std::make_unique<Epetra_Map>(TrilinosWrappers::types::int_type(0),
707  graph = std::make_unique<Epetra_FECrsGraph>(View,
710  0);
711  graph->FillComplete();
712 
713  nonlocal_graph.reset();
714  }
715 
716 
717 
718  void
720  {
721  int ierr;
723  if (nonlocal_graph.get() != nullptr)
724  {
725  if (nonlocal_graph->IndicesAreGlobal() == false &&
726  nonlocal_graph->RowMap().NumMyElements() > 0 &&
728  {
729  // Insert dummy element at (row, column) that corresponds to row 0
730  // in local index counting.
734 
735  // in case we have a square sparsity pattern, add the entry on the
736  // diagonal
739  column = row;
740  // if not, take a column index that we have ourselves since we
741  // know for sure it is there (and it will not create spurious
742  // messages to many ranks like putting index 0 on many processors)
743  else if (column_space_map->NumMyElements() > 0)
745  ierr = nonlocal_graph->InsertGlobalIndices(row, 1, &column);
746  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
747  }
748  Assert(nonlocal_graph->RowMap().NumMyElements() == 0 ||
750  nonlocal_graph->IndicesAreGlobal() == true,
751  ExcInternalError());
752 
753  ierr =
754  nonlocal_graph->FillComplete(*column_space_map, graph->RangeMap());
755  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
756  ierr = nonlocal_graph->OptimizeStorage();
757  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
758  Epetra_Export exporter(nonlocal_graph->RowMap(), graph->RowMap());
759  ierr = graph->Export(*nonlocal_graph, exporter, Add);
760  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
761  ierr = graph->FillComplete(*column_space_map, graph->RangeMap());
762  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
763  }
764  else
765  {
766  // TODO A dynamic_cast fails here, this is suspicious.
767  const auto &range_map =
768  static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
769  ierr = graph->GlobalAssemble(*column_space_map, range_map, true);
770  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
771  }
772 
773  try
774  {
775  ierr = graph->OptimizeStorage();
776  }
777  catch (const int error_code)
778  {
779  AssertThrow(
780  false,
781  ExcMessage(
782  "The Epetra_CrsGraph::OptimizeStorage() function "
783  "has thrown an error with code " +
784  std::to_string(error_code) +
785  ". You will have to look up the exact meaning of this error "
786  "in the Trilinos source code, but oftentimes, this function "
787  "throwing an error indicates that you are trying to allocate "
788  "more than 2,147,483,647 nonzero entries in the sparsity "
789  "pattern on the local process; this will not work because "
790  "Epetra indexes entries with a simple 'signed int'."));
791  }
792  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
793 
794  // Check consistency between the sizes set at the beginning and what
795  // Trilinos stores:
796  using namespace deal_II_exceptions::internals;
798  ExcInternalError());
800  ExcInternalError());
801  }
802 
803 
804 
805  bool
807  {
808  return graph->RowMap().LID(
809  static_cast<TrilinosWrappers::types::int_type>(i)) != -1;
810  }
811 
812 
813 
814  bool
816  {
817  if (!row_is_stored_locally(i))
818  {
819  return false;
820  }
821  else
822  {
823  // Extract local indices in
824  // the matrix.
825  int trilinos_i =
826  graph->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
827  trilinos_j =
828  graph->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
829 
830  // Check whether the matrix
831  // already is transformed to
832  // local indices.
833  if (graph->Filled() == false)
834  {
835  int nnz_present = graph->NumGlobalIndices(i);
836  int nnz_extracted;
838 
839  // Generate the view and make
840  // sure that we have not generated
841  // an error.
842  // TODO: trilinos_i is the local row index -> it is an int but
843  // ExtractGlobalRowView requires trilinos_i to be the global row
844  // index and thus it should be a long long int
845  int ierr = graph->ExtractGlobalRowView(trilinos_i,
846  nnz_extracted,
847  col_indices);
848  (void)ierr;
849  Assert(ierr == 0, ExcTrilinosError(ierr));
850  Assert(nnz_present == nnz_extracted,
851  ExcDimensionMismatch(nnz_present, nnz_extracted));
852 
853  // Search the index
854  const std::ptrdiff_t local_col_index =
855  std::find(col_indices, col_indices + nnz_present, trilinos_j) -
856  col_indices;
857 
858  if (local_col_index == nnz_present)
859  return false;
860  }
861  else
862  {
863  // Prepare pointers for extraction
864  // of a view of the row.
865  int nnz_present = graph->NumGlobalIndices(i);
866  int nnz_extracted;
867  int *col_indices;
868 
869  // Generate the view and make
870  // sure that we have not generated
871  // an error.
872  int ierr =
873  graph->ExtractMyRowView(trilinos_i, nnz_extracted, col_indices);
874  (void)ierr;
875  Assert(ierr == 0, ExcTrilinosError(ierr));
876 
877  Assert(nnz_present == nnz_extracted,
878  ExcDimensionMismatch(nnz_present, nnz_extracted));
879 
880  // Search the index
881  const std::ptrdiff_t local_col_index =
882  std::find(col_indices, col_indices + nnz_present, trilinos_j) -
883  col_indices;
884 
885  if (local_col_index == nnz_present)
886  return false;
887  }
888  }
889 
890  return true;
891  }
892 
893 
894 
897  {
898  size_type local_b = 0;
900  for (int i = 0; i < static_cast<int>(local_size()); ++i)
901  {
902  int *indices;
903  int num_entries;
904  graph->ExtractMyRowView(i, num_entries, indices);
905  for (unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
906  ++j)
907  {
908  if (static_cast<size_type>(std::abs(i - indices[j])) > local_b)
909  local_b = std::abs(i - indices[j]);
910  }
911  }
912  graph->Comm().MaxAll(reinterpret_cast<TrilinosWrappers::types::int_type *>(
913  &local_b),
914  &global_b,
915  1);
916  return static_cast<size_type>(global_b);
917  }
918 
919 
920 
921  unsigned int
923  {
924  int n_rows = graph->NumMyRows();
925 
926  return n_rows;
927  }
928 
929 
930 
931  std::pair<SparsityPattern::size_type, SparsityPattern::size_type>
933  {
936  end = TrilinosWrappers::max_my_gid(graph->RowMap()) + 1;
937 
938  return std::make_pair(begin, end);
939  }
940 
941 
942 
943  std::uint64_t
945  {
947 
948  return static_cast<std::uint64_t>(nnz);
949  }
950 
951 
952 
953  unsigned int
955  {
956  int nnz = graph->MaxNumIndices();
957 
958  return static_cast<unsigned int>(nnz);
959  }
960 
961 
962 
965  {
966  Assert(row < n_rows(), ExcInternalError());
967 
968  // Get a representation of the where the present row is located on
969  // the current processor
971  graph->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
972 
973  // On the processor who owns this row, we'll have a non-negative
974  // value for `local_row` and can ask for the length of the row.
975  if (local_row >= 0)
976  return graph->NumMyIndices(local_row);
977  else
978  return static_cast<size_type>(-1);
979  }
980 
981 
982 
983  void
985  const ArrayView<const size_type> &columns,
986  const bool indices_are_sorted)
987  {
988  add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
989  }
990 
991 
992 
993  const Epetra_Map &
995  {
996  // TODO A dynamic_cast fails here, this is suspicious.
997  const auto &domain_map =
998  static_cast<const Epetra_Map &>(graph->DomainMap()); // NOLINT
999  return domain_map;
1000  }
1001 
1002 
1003 
1004  const Epetra_Map &
1006  {
1007  // TODO A dynamic_cast fails here, this is suspicious.
1008  const auto &range_map =
1009  static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
1010  return range_map;
1011  }
1012 
1013 
1014 
1015  MPI_Comm
1017  {
1018  const Epetra_MpiComm *mpi_comm =
1019  dynamic_cast<const Epetra_MpiComm *>(&graph->RangeMap().Comm());
1020  Assert(mpi_comm != nullptr, ExcInternalError());
1021  return mpi_comm->Comm();
1022  }
1023 
1024 
1025 
1026  void
1028  {
1029  Assert(false, ExcNotImplemented());
1030  }
1031 
1032 
1033 
1034  // As of now, no particularly neat
1035  // output is generated in case of
1036  // multiple processors.
1037  void
1038  SparsityPattern::print(std::ostream &out,
1039  const bool write_extended_trilinos_info) const
1040  {
1041  if (write_extended_trilinos_info)
1042  out << *graph;
1043  else
1044  {
1045  int *indices;
1046  int num_entries;
1047 
1048  for (int i = 0; i < graph->NumMyRows(); ++i)
1049  {
1050  graph->ExtractMyRowView(i, num_entries, indices);
1051  for (int j = 0; j < num_entries; ++j)
1052  out << "(" << TrilinosWrappers::global_index(graph->RowMap(), i)
1053  << ","
1054  << TrilinosWrappers::global_index(graph->ColMap(), indices[j])
1055  << ") " << std::endl;
1056  }
1057  }
1058 
1059  AssertThrow(out.fail() == false, ExcIO());
1060  }
1061 
1062 
1063 
1064  void
1065  SparsityPattern::print_gnuplot(std::ostream &out) const
1066  {
1067  Assert(graph->Filled() == true, ExcInternalError());
1068  for (::types::global_dof_index row = 0; row < local_size(); ++row)
1069  {
1070  int *indices;
1071  int num_entries;
1072  graph->ExtractMyRowView(row, num_entries, indices);
1073 
1074  Assert(num_entries >= 0, ExcInternalError());
1075  // avoid sign comparison warning
1076  const ::types::global_dof_index num_entries_ = num_entries;
1077  for (::types::global_dof_index j = 0; j < num_entries_; ++j)
1078  // while matrix entries are usually
1079  // written (i,j), with i vertical and
1080  // j horizontal, gnuplot output is
1081  // x-y, that is we have to exchange
1082  // the order of output
1083  out << static_cast<int>(
1084  TrilinosWrappers::global_index(graph->ColMap(), indices[j]))
1085  << " "
1086  << -static_cast<int>(
1087  TrilinosWrappers::global_index(graph->RowMap(), row))
1088  << std::endl;
1089  }
1090 
1091  AssertThrow(out.fail() == false, ExcIO());
1092  }
1093 
1094  // TODO: Implement!
1095  std::size_t
1097  {
1098  Assert(false, ExcNotImplemented());
1099  return 0;
1100  }
1101 
1102 
1103 # ifndef DOXYGEN
1104  // explicit instantiations
1105  //
1106  template void
1107  SparsityPattern::copy_from(const ::SparsityPattern &);
1108  template void
1109  SparsityPattern::copy_from(const ::DynamicSparsityPattern &);
1110 
1111  template void
1113  const ::SparsityPattern &,
1114  const MPI_Comm,
1115  bool);
1116  template void
1118  const ::DynamicSparsityPattern &,
1119  const MPI_Comm,
1120  bool);
1121 
1122 
1123  template void
1125  const IndexSet &,
1126  const ::SparsityPattern &,
1127  const MPI_Comm,
1128  bool);
1129  template void
1131  const IndexSet &,
1132  const ::DynamicSparsityPattern &,
1133  const MPI_Comm,
1134  bool);
1135 # endif
1136 
1137 } // namespace TrilinosWrappers
1138 
1140 
1141 #endif // DEAL_II_WITH_TRILINOS
iterator begin() const
Definition: array_view.h:591
iterator end() const
Definition: array_view.h:600
size_type size() const
Definition: index_set.h:1696
size_type n_elements() const
Definition: index_set.h:1854
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:833
void subtract_set(const IndexSet &other)
Definition: index_set.cc:288
virtual void resize(const size_type rows, const size_type cols)
types::global_dof_index size_type
size_type n_rows() const
size_type n_cols() const
std::shared_ptr< const std::vector< size_type > > colnum_cache
size_type row_length(const size_type row) const
std::unique_ptr< Epetra_FECrsGraph > graph
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
void print_gnuplot(std::ostream &out) const
const_iterator end() const
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false) override
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
bool exists(const size_type i, const size_type j) const
const Epetra_Map & domain_partitioner() const
std::pair< size_type, size_type > local_range() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
void copy_from(const SparsityPattern &input_sparsity_pattern)
std::unique_ptr< Epetra_Map > column_space_map
const_iterator begin() const
SparsityPattern & operator=(const SparsityPattern &input_sparsity_pattern)
bool row_is_stored_locally(const size_type i) const
void reinit(const size_type m, const size_type n, const size_type n_entries_per_row=0)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1124
typename ::internal::FEInterfaceViews::ViewType< dim, spacedim, Extractor >::type View
types::global_dof_index size_type
Definition: cuda_kernels.h:45
std::string to_string(const T &t)
Definition: patterns.h:2391
long long int int64_type
Definition: types.h:185
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
TrilinosWrappers::types::int_type n_global_rows(const Epetra_CrsGraph &graph)
TrilinosWrappers::types::int_type min_my_gid(const Epetra_BlockMap &map)
TrilinosWrappers::types::int64_type n_global_entries(const Epetra_CrsGraph &graph)
TrilinosWrappers::types::int64_type n_global_elements(const Epetra_BlockMap &map)
TrilinosWrappers::types::int_type max_my_gid(const Epetra_BlockMap &map)
TrilinosWrappers::types::int_type n_global_cols(const Epetra_CrsGraph &graph)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition: mpi.cc:149
const Epetra_Comm & comm_self()
constexpr bool compare_for_equality(const T &t, const U &u)
Definition: exceptions.h:1746
Definition: types.h:33
unsigned int global_dof_index
Definition: types.h:82