Reference documentation for deal.II version GIT dad323def1 2022-06-25 19: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 - 2022 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  {
86  std::make_unique<Epetra_Map>(TrilinosWrappers::types::int_type(0),
89  graph = std::make_unique<Epetra_FECrsGraph>(View,
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 
261  TrilinosWrappers::min_my_gid(row_map)) *
262  std::uint64_t(n_entries_per_row) <
263  static_cast<std::uint64_t>(std::numeric_limits<int>::max()),
264  ExcMessage("The TrilinosWrappers use Epetra internally which "
265  "uses 'signed int' to represent local indices. "
266  "Therefore, only 2,147,483,647 nonzero matrix "
267  "entries can be stored on a single process, "
268  "but you are requesting more than that. "
269  "If possible, use more MPI processes."));
270 
271 
272  // for more than one processor, need to specify only row map first and
273  // let the matrix entries decide about the column map (which says which
274  // columns are present in the matrix, not to be confused with the
275  // col_map that tells how the domain dofs of the matrix will be
276  // distributed). for only one processor, we can directly assign the
277  // columns as well. If we use a recent Trilinos version, we can also
278  // require building a non-local graph which gives us thread-safe
279  // initialization.
280  if (row_map.Comm().NumProc() > 1)
281  graph = std::make_unique<Epetra_FECrsGraph>(
282  Copy, row_map, n_entries_per_row, false
283  // TODO: Check which new Trilinos version supports this...
284  // Remember to change tests/trilinos/assemble_matrix_parallel_07, too.
285  //#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
286  // , true
287  //#endif
288  );
289  else
290  graph = std::make_unique<Epetra_FECrsGraph>(
291  Copy, row_map, col_map, n_entries_per_row, false);
292  }
293 
294 
295 
296  void
297  reinit_sp(const Epetra_Map & row_map,
298  const Epetra_Map & col_map,
299  const std::vector<size_type> & n_entries_per_row,
300  std::unique_ptr<Epetra_Map> & column_space_map,
301  std::unique_ptr<Epetra_FECrsGraph> &graph,
302  std::unique_ptr<Epetra_CrsGraph> & nonlocal_graph)
303  {
304  Assert(row_map.IsOneToOne(),
305  ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
306  "the maps of different processors."));
307  Assert(col_map.IsOneToOne(),
308  ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
309  "the maps of different processors."));
310 
311  // release memory before reallocation
312  nonlocal_graph.reset();
313  graph.reset();
314  AssertDimension(n_entries_per_row.size(),
316 
317  column_space_map = std::make_unique<Epetra_Map>(col_map);
318 
319  // Translate the vector of row lengths into one that only stores
320  // those entries that related to the locally stored rows of the matrix:
321  std::vector<int> local_entries_per_row(
324  for (unsigned int i = 0; i < local_entries_per_row.size(); ++i)
325  local_entries_per_row[i] =
326  n_entries_per_row[TrilinosWrappers::min_my_gid(row_map) + i];
327 
328  AssertThrow(std::accumulate(local_entries_per_row.begin(),
329  local_entries_per_row.end(),
330  std::uint64_t(0)) <
331  static_cast<std::uint64_t>(std::numeric_limits<int>::max()),
332  ExcMessage("The TrilinosWrappers use Epetra internally which "
333  "uses 'signed int' to represent local indices. "
334  "Therefore, only 2,147,483,647 nonzero matrix "
335  "entries can be stored on a single process, "
336  "but you are requesting more than that. "
337  "If possible, use more MPI processes."));
338 
339  if (row_map.Comm().NumProc() > 1)
340  graph = std::make_unique<Epetra_FECrsGraph>(
341  Copy, row_map, local_entries_per_row.data(), false
342  // TODO: Check which new Trilinos version supports this...
343  // Remember to change tests/trilinos/assemble_matrix_parallel_07, too.
344  //#if DEAL_II_TRILINOS_VERSION_GTE(11,14,0)
345  // , true
346  //#endif
347  );
348  else
349  graph = std::make_unique<Epetra_FECrsGraph>(
350  Copy, row_map, col_map, local_entries_per_row.data(), false);
351  }
352 
353 
354 
355  template <typename SparsityPatternType>
356  void
357  reinit_sp(const Epetra_Map & row_map,
358  const Epetra_Map & col_map,
359  const SparsityPatternType & sp,
360  const bool exchange_data,
361  std::unique_ptr<Epetra_Map> & column_space_map,
362  std::unique_ptr<Epetra_FECrsGraph> &graph,
363  std::unique_ptr<Epetra_CrsGraph> & nonlocal_graph)
364  {
365  nonlocal_graph.reset();
366  graph.reset();
367 
368  AssertDimension(sp.n_rows(),
370  AssertDimension(sp.n_cols(),
372 
373  column_space_map = std::make_unique<Epetra_Map>(col_map);
374 
375  Assert(row_map.LinearMap() == true,
376  ExcMessage(
377  "This function only works if the row map is contiguous."));
378 
379  const size_type first_row = TrilinosWrappers::min_my_gid(row_map),
380  last_row = TrilinosWrappers::max_my_gid(row_map) + 1;
381  std::vector<int> n_entries_per_row(last_row - first_row);
382 
383  // Trilinos wants the row length as an int this is hopefully never going
384  // to be a problem.
385  for (size_type row = first_row; row < last_row; ++row)
386  n_entries_per_row[row - first_row] =
387  static_cast<int>(sp.row_length(row));
388 
389  AssertThrow(std::accumulate(n_entries_per_row.begin(),
390  n_entries_per_row.end(),
391  std::uint64_t(0)) <
392  static_cast<std::uint64_t>(std::numeric_limits<int>::max()),
393  ExcMessage("The TrilinosWrappers use Epetra internally which "
394  "uses 'signed int' to represent local indices. "
395  "Therefore, only 2,147,483,647 nonzero matrix "
396  "entries can be stored on a single process, "
397  "but you are requesting more than that. "
398  "If possible, use more MPI processes."));
399 
400  if (row_map.Comm().NumProc() > 1)
401  graph = std::make_unique<Epetra_FECrsGraph>(Copy,
402  row_map,
403  n_entries_per_row.data(),
404  false);
405  else
406  graph = std::make_unique<Epetra_FECrsGraph>(
407  Copy, row_map, col_map, n_entries_per_row.data(), false);
408 
409  AssertDimension(sp.n_rows(), n_global_rows(*graph));
410 
411  std::vector<TrilinosWrappers::types::int_type> row_indices;
412 
413  // Include possibility to exchange data since DynamicSparsityPattern is
414  // able to do so
415  if (exchange_data == false)
416  for (size_type row = first_row; row < last_row; ++row)
417  {
418  const TrilinosWrappers::types::int_type row_length =
419  sp.row_length(row);
420  if (row_length == 0)
421  continue;
422 
423  row_indices.resize(row_length, -1);
424  {
425  typename SparsityPatternType::iterator p = sp.begin(row);
426  // avoid incrementing p over the end of the current row because
427  // it is slow for DynamicSparsityPattern in parallel
428  for (int col = 0; col < row_length;)
429  {
430  row_indices[col++] = p->column();
431  if (col < row_length)
432  ++p;
433  }
434  }
435  graph->Epetra_CrsGraph::InsertGlobalIndices(row,
436  row_length,
437  row_indices.data());
438  }
439  else
440  for (size_type row = 0; row < sp.n_rows(); ++row)
441  {
442  const TrilinosWrappers::types::int_type row_length =
443  sp.row_length(row);
444  if (row_length == 0)
445  continue;
446 
447  row_indices.resize(row_length, -1);
448  {
449  typename SparsityPatternType::iterator p = sp.begin(row);
450  // avoid incrementing p over the end of the current row because
451  // it is slow for DynamicSparsityPattern in parallel
452  for (int col = 0; col < row_length;)
453  {
454  row_indices[col++] = p->column();
455  if (col < row_length)
456  ++p;
457  }
458  }
459  const TrilinosWrappers::types::int_type trilinos_row = row;
460  graph->InsertGlobalIndices(1,
461  &trilinos_row,
462  row_length,
463  row_indices.data());
464  }
465 
466  // TODO A dynamic_cast fails here, this is suspicious.
467  const auto &range_map =
468  static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
469  int ierr = graph->GlobalAssemble(*column_space_map, range_map, true);
470  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
471 
472  ierr = graph->OptimizeStorage();
473  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
474  }
475  } // namespace
476 
477 
478 
479  void
480  SparsityPattern::reinit(const IndexSet &parallel_partitioning,
481  const MPI_Comm &communicator,
482  const size_type n_entries_per_row)
483  {
484  Epetra_Map map =
485  parallel_partitioning.make_trilinos_map(communicator, false);
486  reinit_sp(
487  map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
488  }
489 
490 
491 
492  void
493  SparsityPattern::reinit(const IndexSet & parallel_partitioning,
494  const MPI_Comm & communicator,
495  const std::vector<size_type> &n_entries_per_row)
496  {
497  Epetra_Map map =
498  parallel_partitioning.make_trilinos_map(communicator, false);
499  reinit_sp(
500  map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
501  }
502 
503 
504 
505  void
506  SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
507  const IndexSet &col_parallel_partitioning,
508  const MPI_Comm &communicator,
509  const size_type n_entries_per_row)
510  {
511  Epetra_Map row_map =
512  row_parallel_partitioning.make_trilinos_map(communicator, false);
513  Epetra_Map col_map =
514  col_parallel_partitioning.make_trilinos_map(communicator, false);
515  reinit_sp(row_map,
516  col_map,
517  n_entries_per_row,
519  graph,
521  }
522 
523 
524 
525  void
526  SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
527  const IndexSet &col_parallel_partitioning,
528  const MPI_Comm &communicator,
529  const std::vector<size_type> &n_entries_per_row)
530  {
531  Epetra_Map row_map =
532  row_parallel_partitioning.make_trilinos_map(communicator, false);
533  Epetra_Map col_map =
534  col_parallel_partitioning.make_trilinos_map(communicator, false);
535  reinit_sp(row_map,
536  col_map,
537  n_entries_per_row,
539  graph,
541  }
542 
543 
544 
545  void
546  SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
547  const IndexSet &col_parallel_partitioning,
548  const IndexSet &writable_rows,
549  const MPI_Comm &communicator,
550  const size_type n_entries_per_row)
551  {
552  Epetra_Map row_map =
553  row_parallel_partitioning.make_trilinos_map(communicator, false);
554  Epetra_Map col_map =
555  col_parallel_partitioning.make_trilinos_map(communicator, false);
556  reinit_sp(row_map,
557  col_map,
558  n_entries_per_row,
560  graph,
562 
563  IndexSet nonlocal_partitioner = writable_rows;
564  AssertDimension(nonlocal_partitioner.size(),
565  row_parallel_partitioning.size());
566 # ifdef DEBUG
567  {
568  IndexSet tmp = writable_rows & row_parallel_partitioning;
569  Assert(tmp == row_parallel_partitioning,
570  ExcMessage(
571  "The set of writable rows passed to this method does not "
572  "contain the locally owned rows, which is not allowed."));
573  }
574 # endif
575  nonlocal_partitioner.subtract_set(row_parallel_partitioning);
576  if (Utilities::MPI::n_mpi_processes(communicator) > 1)
577  {
578  Epetra_Map nonlocal_map =
579  nonlocal_partitioner.make_trilinos_map(communicator, true);
581  std::make_unique<Epetra_CrsGraph>(Copy, nonlocal_map, 0);
582  }
583  else
584  Assert(nonlocal_partitioner.n_elements() == 0, ExcInternalError());
585  }
586 
587 
588 
589  template <typename SparsityPatternType>
590  void
592  const IndexSet & row_parallel_partitioning,
593  const IndexSet & col_parallel_partitioning,
594  const SparsityPatternType &nontrilinos_sparsity_pattern,
595  const MPI_Comm & communicator,
596  const bool exchange_data)
597  {
598  Epetra_Map row_map =
599  row_parallel_partitioning.make_trilinos_map(communicator, false);
600  Epetra_Map col_map =
601  col_parallel_partitioning.make_trilinos_map(communicator, false);
602  reinit_sp(row_map,
603  col_map,
604  nontrilinos_sparsity_pattern,
605  exchange_data,
607  graph,
609  }
610 
611 
612 
613  template <typename SparsityPatternType>
614  void
616  const IndexSet & parallel_partitioning,
617  const SparsityPatternType &nontrilinos_sparsity_pattern,
618  const MPI_Comm & communicator,
619  const bool exchange_data)
620  {
621  Epetra_Map map =
622  parallel_partitioning.make_trilinos_map(communicator, false);
623  reinit_sp(map,
624  map,
625  nontrilinos_sparsity_pattern,
626  exchange_data,
628  graph,
630  }
631 
632 
633 
636  {
637  Assert(false, ExcNotImplemented());
638  return *this;
639  }
640 
641 
642 
643  void
645  {
646  column_space_map = std::make_unique<Epetra_Map>(*sp.column_space_map);
647  graph = std::make_unique<Epetra_FECrsGraph>(*sp.graph);
648 
649  if (sp.nonlocal_graph.get() != nullptr)
650  nonlocal_graph = std::make_unique<Epetra_CrsGraph>(*sp.nonlocal_graph);
651  else
652  nonlocal_graph.reset();
653  }
654 
655 
656 
657  template <typename SparsityPatternType>
658  void
659  SparsityPattern::copy_from(const SparsityPatternType &sp)
660  {
661  const Epetra_Map rows(TrilinosWrappers::types::int_type(sp.n_rows()),
662  0,
664  const Epetra_Map columns(TrilinosWrappers::types::int_type(sp.n_cols()),
665  0,
667 
668  reinit_sp(
669  rows, columns, sp, false, column_space_map, graph, nonlocal_graph);
670  }
671 
672 
673 
674  void
676  {
677  // When we clear the matrix, reset
678  // the pointer and generate an
679  // empty sparsity pattern.
681  std::make_unique<Epetra_Map>(TrilinosWrappers::types::int_type(0),
684  graph = std::make_unique<Epetra_FECrsGraph>(View,
687  0);
688  graph->FillComplete();
689 
690  nonlocal_graph.reset();
691  }
692 
693 
694 
695  void
697  {
698  int ierr;
700  if (nonlocal_graph.get() != nullptr)
701  {
702  if (nonlocal_graph->IndicesAreGlobal() == false &&
703  nonlocal_graph->RowMap().NumMyElements() > 0 &&
705  {
706  // Insert dummy element at (row, column) that corresponds to row 0
707  // in local index counting.
711 
712  // in case we have a square sparsity pattern, add the entry on the
713  // diagonal
716  column = row;
717  // if not, take a column index that we have ourselves since we
718  // know for sure it is there (and it will not create spurious
719  // messages to many ranks like putting index 0 on many processors)
720  else if (column_space_map->NumMyElements() > 0)
722  ierr = nonlocal_graph->InsertGlobalIndices(row, 1, &column);
723  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
724  }
725  Assert(nonlocal_graph->RowMap().NumMyElements() == 0 ||
727  nonlocal_graph->IndicesAreGlobal() == true,
728  ExcInternalError());
729 
730  ierr =
731  nonlocal_graph->FillComplete(*column_space_map, graph->RangeMap());
732  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
733  ierr = nonlocal_graph->OptimizeStorage();
734  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
735  Epetra_Export exporter(nonlocal_graph->RowMap(), graph->RowMap());
736  ierr = graph->Export(*nonlocal_graph, exporter, Add);
737  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
738  ierr = graph->FillComplete(*column_space_map, graph->RangeMap());
739  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
740  }
741  else
742  {
743  // TODO A dynamic_cast fails here, this is suspicious.
744  const auto &range_map =
745  static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
746  ierr = graph->GlobalAssemble(*column_space_map, range_map, true);
747  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
748  }
749 
750  try
751  {
752  ierr = graph->OptimizeStorage();
753  }
754  catch (const int error_code)
755  {
756  AssertThrow(
757  false,
758  ExcMessage(
759  "The Epetra_CrsGraph::OptimizeStorage() function "
760  "has thrown an error with code " +
761  std::to_string(error_code) +
762  ". You will have to look up the exact meaning of this error "
763  "in the Trilinos source code, but oftentimes, this function "
764  "throwing an error indicates that you are trying to allocate "
765  "more than 2,147,483,647 nonzero entries in the sparsity "
766  "pattern on the local process; this will not work because "
767  "Epetra indexes entries with a simple 'signed int'."));
768  }
769  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
770  }
771 
772 
773 
774  bool
776  {
777  return graph->RowMap().LID(
778  static_cast<TrilinosWrappers::types::int_type>(i)) != -1;
779  }
780 
781 
782 
783  bool
785  {
786  if (!row_is_stored_locally(i))
787  {
788  return false;
789  }
790  else
791  {
792  // Extract local indices in
793  // the matrix.
794  int trilinos_i =
795  graph->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
796  trilinos_j =
797  graph->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
798 
799  // Check whether the matrix
800  // already is transformed to
801  // local indices.
802  if (graph->Filled() == false)
803  {
804  int nnz_present = graph->NumGlobalIndices(i);
805  int nnz_extracted;
807 
808  // Generate the view and make
809  // sure that we have not generated
810  // an error.
811  // TODO: trilinos_i is the local row index -> it is an int but
812  // ExtractGlobalRowView requires trilinos_i to be the global row
813  // index and thus it should be a long long int
814  int ierr = graph->ExtractGlobalRowView(trilinos_i,
815  nnz_extracted,
816  col_indices);
817  (void)ierr;
818  Assert(ierr == 0, ExcTrilinosError(ierr));
819  Assert(nnz_present == nnz_extracted,
820  ExcDimensionMismatch(nnz_present, nnz_extracted));
821 
822  // Search the index
823  const std::ptrdiff_t local_col_index =
824  std::find(col_indices, col_indices + nnz_present, trilinos_j) -
825  col_indices;
826 
827  if (local_col_index == nnz_present)
828  return false;
829  }
830  else
831  {
832  // Prepare pointers for extraction
833  // of a view of the row.
834  int nnz_present = graph->NumGlobalIndices(i);
835  int nnz_extracted;
836  int *col_indices;
837 
838  // Generate the view and make
839  // sure that we have not generated
840  // an error.
841  int ierr =
842  graph->ExtractMyRowView(trilinos_i, nnz_extracted, col_indices);
843  (void)ierr;
844  Assert(ierr == 0, ExcTrilinosError(ierr));
845 
846  Assert(nnz_present == nnz_extracted,
847  ExcDimensionMismatch(nnz_present, nnz_extracted));
848 
849  // Search the index
850  const std::ptrdiff_t local_col_index =
851  std::find(col_indices, col_indices + nnz_present, trilinos_j) -
852  col_indices;
853 
854  if (local_col_index == nnz_present)
855  return false;
856  }
857  }
858 
859  return true;
860  }
861 
862 
863 
866  {
867  size_type local_b = 0;
869  for (int i = 0; i < static_cast<int>(local_size()); ++i)
870  {
871  int *indices;
872  int num_entries;
873  graph->ExtractMyRowView(i, num_entries, indices);
874  for (unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
875  ++j)
876  {
877  if (static_cast<size_type>(std::abs(i - indices[j])) > local_b)
878  local_b = std::abs(i - indices[j]);
879  }
880  }
881  graph->Comm().MaxAll(reinterpret_cast<TrilinosWrappers::types::int_type *>(
882  &local_b),
883  &global_b,
884  1);
885  return static_cast<size_type>(global_b);
886  }
887 
888 
889 
892  {
894  return n_rows;
895  }
896 
897 
898 
901  {
903  if (graph->Filled() == true)
905  else
907 
908  return n_cols;
909  }
910 
911 
912 
913  unsigned int
915  {
916  int n_rows = graph->NumMyRows();
917 
918  return n_rows;
919  }
920 
921 
922 
923  std::pair<SparsityPattern::size_type, SparsityPattern::size_type>
925  {
928  end = TrilinosWrappers::max_my_gid(graph->RowMap()) + 1;
929 
930  return std::make_pair(begin, end);
931  }
932 
933 
934 
935  std::uint64_t
937  {
939 
940  return static_cast<std::uint64_t>(nnz);
941  }
942 
943 
944 
945  unsigned int
947  {
948  int nnz = graph->MaxNumIndices();
949 
950  return static_cast<unsigned int>(nnz);
951  }
952 
953 
954 
957  {
958  Assert(row < n_rows(), ExcInternalError());
959 
960  // Get a representation of the where the present row is located on
961  // the current processor
963  graph->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
964 
965  // On the processor who owns this row, we'll have a non-negative
966  // value for `local_row` and can ask for the length of the row.
967  if (local_row >= 0)
968  return graph->NumMyIndices(local_row);
969  else
970  return static_cast<size_type>(-1);
971  }
972 
973 
974 
975  const Epetra_Map &
977  {
978  // TODO A dynamic_cast fails here, this is suspicious.
979  const auto &domain_map =
980  static_cast<const Epetra_Map &>(graph->DomainMap()); // NOLINT
981  return domain_map;
982  }
983 
984 
985 
986  const Epetra_Map &
988  {
989  // TODO A dynamic_cast fails here, this is suspicious.
990  const auto &range_map =
991  static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
992  return range_map;
993  }
994 
995 
996 
997  MPI_Comm
999  {
1000  const Epetra_MpiComm *mpi_comm =
1001  dynamic_cast<const Epetra_MpiComm *>(&graph->RangeMap().Comm());
1002  Assert(mpi_comm != nullptr, ExcInternalError());
1003  return mpi_comm->Comm();
1004  }
1005 
1006 
1007 
1008  void
1010  {
1011  Assert(false, ExcNotImplemented());
1012  }
1013 
1014 
1015 
1016  // As of now, no particularly neat
1017  // output is generated in case of
1018  // multiple processors.
1019  void
1020  SparsityPattern::print(std::ostream &out,
1021  const bool write_extended_trilinos_info) const
1022  {
1023  if (write_extended_trilinos_info)
1024  out << *graph;
1025  else
1026  {
1027  int *indices;
1028  int num_entries;
1029 
1030  for (int i = 0; i < graph->NumMyRows(); ++i)
1031  {
1032  graph->ExtractMyRowView(i, num_entries, indices);
1033  for (int j = 0; j < num_entries; ++j)
1034  out << "(" << TrilinosWrappers::global_index(graph->RowMap(), i)
1035  << ","
1036  << TrilinosWrappers::global_index(graph->ColMap(), indices[j])
1037  << ") " << std::endl;
1038  }
1039  }
1040 
1041  AssertThrow(out.fail() == false, ExcIO());
1042  }
1043 
1044 
1045 
1046  void
1047  SparsityPattern::print_gnuplot(std::ostream &out) const
1048  {
1049  Assert(graph->Filled() == true, ExcInternalError());
1050  for (::types::global_dof_index row = 0; row < local_size(); ++row)
1051  {
1052  int *indices;
1053  int num_entries;
1054  graph->ExtractMyRowView(row, num_entries, indices);
1055 
1056  Assert(num_entries >= 0, ExcInternalError());
1057  // avoid sign comparison warning
1058  const ::types::global_dof_index num_entries_ = num_entries;
1059  for (::types::global_dof_index j = 0; j < num_entries_; ++j)
1060  // while matrix entries are usually
1061  // written (i,j), with i vertical and
1062  // j horizontal, gnuplot output is
1063  // x-y, that is we have to exchange
1064  // the order of output
1065  out << static_cast<int>(
1066  TrilinosWrappers::global_index(graph->ColMap(), indices[j]))
1067  << " "
1068  << -static_cast<int>(
1069  TrilinosWrappers::global_index(graph->RowMap(), row))
1070  << std::endl;
1071  }
1072 
1073  AssertThrow(out.fail() == false, ExcIO());
1074  }
1075 
1076  // TODO: Implement!
1077  std::size_t
1079  {
1080  Assert(false, ExcNotImplemented());
1081  return 0;
1082  }
1083 
1084 
1085 # ifndef DOXYGEN
1086  // explicit instantiations
1087  //
1088  template void
1089  SparsityPattern::copy_from(const ::SparsityPattern &);
1090  template void
1091  SparsityPattern::copy_from(const ::DynamicSparsityPattern &);
1092 
1093  template void
1095  const ::SparsityPattern &,
1096  const MPI_Comm &,
1097  bool);
1098  template void
1100  const ::DynamicSparsityPattern &,
1101  const MPI_Comm &,
1102  bool);
1103 
1104 
1105  template void
1107  const IndexSet &,
1108  const ::SparsityPattern &,
1109  const MPI_Comm &,
1110  bool);
1111  template void
1113  const IndexSet &,
1114  const ::DynamicSparsityPattern &,
1115  const MPI_Comm &,
1116  bool);
1117 # endif
1118 
1119 } // namespace TrilinosWrappers
1120 
1122 
1123 #endif // DEAL_II_WITH_TRILINOS
size_type size() const
Definition: index_set.h:1635
size_type n_elements() const
Definition: index_set.h:1783
void subtract_set(const IndexSet &other)
Definition: index_set.cc:266
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:767
std::shared_ptr< const std::vector< size_type > > colnum_cache
size_type row_length(const size_type row) const
const_iterator end() const
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
void copy_from(const SparsityPattern &input_sparsity_pattern)
const Epetra_Map & range_partitioner() const
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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcInternalError()
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
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::string to_string(const T &t)
Definition: patterns.h:2403
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
std::unique_ptr< Epetra_Map > column_space_map
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1063
typename ::internal::FEInterfaceViews::ViewType< dim, spacedim, Extractor >::type View
types::global_dof_index size_type
Definition: cuda_kernels.h:45
long long int int64_type
Definition: types.h:173
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:145
const Epetra_Comm & comm_self()
Definition: utilities.cc:1110
Definition: types.h:32
unsigned int global_dof_index
Definition: types.h:76