Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.4.1
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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>
23
26
27# include <Epetra_Export.h>
28
30
31namespace TrilinosWrappers
32{
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),
131 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,
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 {
239 using size_type = SparsityPattern::size_type;
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
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,
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,
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,
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 {
757 false,
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
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:1636
size_type n_elements() const
Definition: index_set.h:1834
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:601
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 & ExcIO()
std::unique_ptr< Epetra_FECrsGraph > graph
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
static ::ExceptionBase & ExcNotImplemented()
void print_gnuplot(std::ostream &out) const
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
std::unique_ptr< Epetra_Map > column_space_map
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
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:1014
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:140
const Epetra_Comm & comm_self()
Definition: utilities.cc:1110
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: types.h:32