Reference documentation for deal.II version 9.6.0
\(\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\}}\)
Loading...
Searching...
No Matches
trilinos_sparsity_pattern.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17
18#ifdef DEAL_II_WITH_TRILINOS
19
20# include <deal.II/base/mpi.h>
22
25
26# include <Epetra_Export.h>
27
28# include <limits>
29
31
32namespace TrilinosWrappers
33{
35 {
36 void
38 {
39 // if we are asked to visit the past-the-end line, then simply
40 // release all our caches and go on with life
41 if (this->a_row == sparsity_pattern->n_rows() ||
42 (sparsity_pattern->in_local_range(this->a_row) == false))
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),
133 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,
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 {
241 using size_type = SparsityPattern::size_type;
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
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,
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 for (size_type row = first_row; row < last_row; ++row)
416 {
417 const TrilinosWrappers::types::int_type row_length =
418 sp.row_length(row);
419 if (row_length == 0)
420 continue;
421
422 row_indices.resize(row_length, -1);
423 {
424 typename SparsityPatternType::iterator p = sp.begin(row);
425 // avoid incrementing p over the end of the current row because
426 // it is slow for DynamicSparsityPattern in parallel
427 for (int col = 0; col < row_length;)
428 {
429 row_indices[col++] = p->column();
430 if (col < row_length)
431 ++p;
432 }
433 }
434 if (exchange_data == false)
435 graph->Epetra_CrsGraph::InsertGlobalIndices(row,
436 row_length,
437 row_indices.data());
438 else
439 {
440 // Include possibility to exchange data since
441 // DynamicSparsityPattern is able to do so
442 const TrilinosWrappers::types::int_type trilinos_row = row;
443 graph->InsertGlobalIndices(1,
444 &trilinos_row,
445 row_length,
446 row_indices.data());
447 }
448 }
449
450 // TODO A dynamic_cast fails here, this is suspicious.
451 const auto &range_map =
452 static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
453 int ierr = graph->GlobalAssemble(*column_space_map, range_map, true);
454 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
455
456 ierr = graph->OptimizeStorage();
457 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
458 }
459 } // namespace
460
461
462
463 void
464 SparsityPattern::reinit(const IndexSet &parallel_partitioning,
465 const MPI_Comm communicator,
466 const size_type n_entries_per_row)
467 {
468 SparsityPatternBase::resize(parallel_partitioning.size(),
469 parallel_partitioning.size());
470 Epetra_Map map =
471 parallel_partitioning.make_trilinos_map(communicator, false);
472 reinit_sp(
473 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
474 }
475
476
477
478 void
479 SparsityPattern::reinit(const IndexSet &parallel_partitioning,
480 const MPI_Comm communicator,
481 const std::vector<size_type> &n_entries_per_row)
482 {
483 SparsityPatternBase::resize(parallel_partitioning.size(),
484 parallel_partitioning.size());
485 Epetra_Map map =
486 parallel_partitioning.make_trilinos_map(communicator, false);
487 reinit_sp(
488 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
489 }
490
491
492
493 void
494 SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
495 const IndexSet &col_parallel_partitioning,
496 const MPI_Comm communicator,
497 const size_type n_entries_per_row)
498 {
499 SparsityPatternBase::resize(row_parallel_partitioning.size(),
500 col_parallel_partitioning.size());
501 Epetra_Map row_map =
502 row_parallel_partitioning.make_trilinos_map(communicator, false);
503 Epetra_Map col_map =
504 col_parallel_partitioning.make_trilinos_map(communicator, false);
505 reinit_sp(row_map,
506 col_map,
507 n_entries_per_row,
509 graph,
511 }
512
513
514
515 void
516 SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
517 const IndexSet &col_parallel_partitioning,
518 const MPI_Comm communicator,
519 const std::vector<size_type> &n_entries_per_row)
520 {
521 SparsityPatternBase::resize(row_parallel_partitioning.size(),
522 col_parallel_partitioning.size());
523 Epetra_Map row_map =
524 row_parallel_partitioning.make_trilinos_map(communicator, false);
525 Epetra_Map col_map =
526 col_parallel_partitioning.make_trilinos_map(communicator, false);
527 reinit_sp(row_map,
528 col_map,
529 n_entries_per_row,
531 graph,
533 }
534
535
536
537 void
538 SparsityPattern::reinit(const IndexSet &row_parallel_partitioning,
539 const IndexSet &col_parallel_partitioning,
540 const IndexSet &writable_rows,
541 const MPI_Comm communicator,
542 const size_type n_entries_per_row)
543 {
544 SparsityPatternBase::resize(row_parallel_partitioning.size(),
545 col_parallel_partitioning.size());
546 Epetra_Map row_map =
547 row_parallel_partitioning.make_trilinos_map(communicator, false);
548 Epetra_Map col_map =
549 col_parallel_partitioning.make_trilinos_map(communicator, false);
550 reinit_sp(row_map,
551 col_map,
552 n_entries_per_row,
554 graph,
556
557 IndexSet nonlocal_partitioner = writable_rows;
558 AssertDimension(nonlocal_partitioner.size(),
559 row_parallel_partitioning.size());
560# ifdef DEBUG
561 {
562 IndexSet tmp = writable_rows & row_parallel_partitioning;
563 Assert(tmp == row_parallel_partitioning,
565 "The set of writable rows passed to this method does not "
566 "contain the locally owned rows, which is not allowed."));
567 }
568# endif
569 nonlocal_partitioner.subtract_set(row_parallel_partitioning);
570 if (Utilities::MPI::n_mpi_processes(communicator) > 1)
571 {
572 Epetra_Map nonlocal_map =
573 nonlocal_partitioner.make_trilinos_map(communicator, true);
575 std::make_unique<Epetra_CrsGraph>(Copy, nonlocal_map, 0);
576 }
577 else
578 Assert(nonlocal_partitioner.n_elements() == 0, ExcInternalError());
579 }
580
581
582
583 template <typename SparsityPatternType>
584 void
586 const IndexSet &row_parallel_partitioning,
587 const IndexSet &col_parallel_partitioning,
588 const SparsityPatternType &nontrilinos_sparsity_pattern,
589 const MPI_Comm communicator,
590 const bool exchange_data)
591 {
592 SparsityPatternBase::resize(row_parallel_partitioning.size(),
593 col_parallel_partitioning.size());
594 Epetra_Map row_map =
595 row_parallel_partitioning.make_trilinos_map(communicator, false);
596 Epetra_Map col_map =
597 col_parallel_partitioning.make_trilinos_map(communicator, false);
598 reinit_sp(row_map,
599 col_map,
600 nontrilinos_sparsity_pattern,
601 exchange_data,
603 graph,
605 }
606
607
608
609 template <typename SparsityPatternType>
610 void
612 const IndexSet &parallel_partitioning,
613 const SparsityPatternType &nontrilinos_sparsity_pattern,
614 const MPI_Comm communicator,
615 const bool exchange_data)
616 {
617 AssertDimension(nontrilinos_sparsity_pattern.n_rows(),
618 parallel_partitioning.size());
619 AssertDimension(nontrilinos_sparsity_pattern.n_cols(),
620 parallel_partitioning.size());
621 SparsityPatternBase::resize(parallel_partitioning.size(),
622 parallel_partitioning.size());
623 Epetra_Map map =
624 parallel_partitioning.make_trilinos_map(communicator, false);
625 reinit_sp(map,
626 map,
627 nontrilinos_sparsity_pattern,
628 exchange_data,
630 graph,
632 }
633
634
635
638 {
640 return *this;
641 }
642
643
644
645 void
647 {
649 column_space_map = std::make_unique<Epetra_Map>(*sp.column_space_map);
650 graph = std::make_unique<Epetra_FECrsGraph>(*sp.graph);
651
652 if (sp.nonlocal_graph.get() != nullptr)
653 nonlocal_graph = std::make_unique<Epetra_CrsGraph>(*sp.nonlocal_graph);
654 else
655 nonlocal_graph.reset();
656 }
657
658
659
660 template <typename SparsityPatternType>
661 void
662 SparsityPattern::copy_from(const SparsityPatternType &sp)
663 {
664 SparsityPatternBase::resize(sp.n_rows(), sp.n_cols());
665 const Epetra_Map rows(TrilinosWrappers::types::int_type(sp.n_rows()),
666 0,
668 const Epetra_Map columns(TrilinosWrappers::types::int_type(sp.n_cols()),
669 0,
671
672 reinit_sp(
673 rows, columns, sp, false, column_space_map, graph, nonlocal_graph);
674 }
675
676
677
678 void
680 {
682 // When we clear the matrix, reset
683 // the pointer and generate an
684 // empty sparsity pattern.
686 std::make_unique<Epetra_Map>(TrilinosWrappers::types::int_type(0),
689 graph = std::make_unique<Epetra_FECrsGraph>(View,
692 0);
693 graph->FillComplete();
694
695 nonlocal_graph.reset();
696 }
697
698
699
700 void
702 {
703 int ierr;
705 if (nonlocal_graph.get() != nullptr)
706 {
707 if (nonlocal_graph->IndicesAreGlobal() == false &&
708 nonlocal_graph->RowMap().NumMyElements() > 0 &&
710 {
711 // Insert dummy element at (row, column) that corresponds to row 0
712 // in local index counting.
716
717 // in case we have a square sparsity pattern, add the entry on the
718 // diagonal
721 column = row;
722 // if not, take a column index that we have ourselves since we
723 // know for sure it is there (and it will not create spurious
724 // messages to many ranks like putting index 0 on many processors)
725 else if (column_space_map->NumMyElements() > 0)
727 ierr = nonlocal_graph->InsertGlobalIndices(row, 1, &column);
728 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
729 }
730 Assert(nonlocal_graph->RowMap().NumMyElements() == 0 ||
732 nonlocal_graph->IndicesAreGlobal() == true,
734
735 ierr =
736 nonlocal_graph->FillComplete(*column_space_map, graph->RangeMap());
737 AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
738 ierr = nonlocal_graph->OptimizeStorage();
739 AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
740 Epetra_Export exporter(nonlocal_graph->RowMap(), graph->RowMap());
741 ierr = graph->Export(*nonlocal_graph, exporter, Add);
742 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
743 ierr = graph->FillComplete(*column_space_map, graph->RangeMap());
744 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
745 }
746 else
747 {
748 // TODO A dynamic_cast fails here, this is suspicious.
749 const auto &range_map =
750 static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
751 ierr = graph->GlobalAssemble(*column_space_map, range_map, true);
752 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
753 }
754
755 try
756 {
757 ierr = graph->OptimizeStorage();
758 }
759 catch (const int error_code)
760 {
762 false,
764 "The Epetra_CrsGraph::OptimizeStorage() function "
765 "has thrown an error with code " +
766 std::to_string(error_code) +
767 ". You will have to look up the exact meaning of this error "
768 "in the Trilinos source code, but oftentimes, this function "
769 "throwing an error indicates that you are trying to allocate "
770 "more than 2,147,483,647 nonzero entries in the sparsity "
771 "pattern on the local process; this will not work because "
772 "Epetra indexes entries with a simple 'signed int'."));
773 }
774 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
775
776 // Check consistency between the sizes set at the beginning and what
777 // Trilinos stores:
778 using namespace deal_II_exceptions::internals;
779 Assert(compare_for_equality(n_rows(), n_global_rows(*graph)),
781 Assert(compare_for_equality(n_cols(), n_global_cols(*graph)),
783 }
784
785
786
787 bool
789 {
790 return graph->RowMap().LID(
791 static_cast<TrilinosWrappers::types::int_type>(i)) != -1;
792 }
793
794
795
796 bool
798 {
799 if (!row_is_stored_locally(i))
800 {
801 return false;
802 }
803 else
804 {
805 // Extract local indices in
806 // the matrix.
807 int trilinos_i =
808 graph->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
809 trilinos_j =
810 graph->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
811
812 // Check whether the matrix
813 // already is transformed to
814 // local indices.
815 if (graph->Filled() == false)
816 {
817 int nnz_present = graph->NumGlobalIndices(i);
818 int nnz_extracted;
820
821 // Generate the view and make
822 // sure that we have not generated
823 // an error.
824 // TODO: trilinos_i is the local row index -> it is an int but
825 // ExtractGlobalRowView requires trilinos_i to be the global row
826 // index and thus it should be a long long int
827 int ierr = graph->ExtractGlobalRowView(trilinos_i,
828 nnz_extracted,
829 col_indices);
830 (void)ierr;
831 Assert(ierr == 0, ExcTrilinosError(ierr));
832 Assert(nnz_present == nnz_extracted,
833 ExcDimensionMismatch(nnz_present, nnz_extracted));
834
835 // Search the index
836 const std::ptrdiff_t local_col_index =
837 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
838 col_indices;
839
840 if (local_col_index == nnz_present)
841 return false;
842 }
843 else
844 {
845 // Prepare pointers for extraction
846 // of a view of the row.
847 int nnz_present = graph->NumGlobalIndices(i);
848 int nnz_extracted;
849 int *col_indices;
850
851 // Generate the view and make
852 // sure that we have not generated
853 // an error.
854 int ierr =
855 graph->ExtractMyRowView(trilinos_i, nnz_extracted, col_indices);
856 (void)ierr;
857 Assert(ierr == 0, ExcTrilinosError(ierr));
858
859 Assert(nnz_present == nnz_extracted,
860 ExcDimensionMismatch(nnz_present, nnz_extracted));
861
862 // Search the index
863 const std::ptrdiff_t local_col_index =
864 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
865 col_indices;
866
867 if (local_col_index == nnz_present)
868 return false;
869 }
870 }
871
872 return true;
873 }
874
875
876
879 {
880 size_type local_b = 0;
881 for (int i = 0; i < static_cast<int>(local_size()); ++i)
882 {
883 int *indices;
884 int num_entries;
885 graph->ExtractMyRowView(i, num_entries, indices);
886 for (unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
887 ++j)
888 {
889 if (static_cast<size_type>(std::abs(i - indices[j])) > local_b)
890 local_b = std::abs(i - indices[j]);
891 }
892 }
893
895 graph->Comm().MaxAll(reinterpret_cast<TrilinosWrappers::types::int_type *>(
896 &local_b),
897 &global_b,
898 1);
899 return static_cast<size_type>(global_b);
900 }
901
902
903
904 unsigned int
906 {
907 return graph->NumMyRows();
908 }
909
910
911
912 std::pair<SparsityPattern::size_type, SparsityPattern::size_type>
914 {
916 const size_type end = TrilinosWrappers::max_my_gid(graph->RowMap()) + 1;
917
918 return {begin, end};
919 }
920
921
922
923 std::uint64_t
928
929
930
931 unsigned int
933 {
934 return graph->MaxNumIndices();
935 }
936
937
938
941 {
942 Assert(row < n_rows(), ExcInternalError());
943
944 // Get a representation of the where the present row is located on
945 // the current processor
947 graph->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
948
949 // On the processor who owns this row, we'll have a non-negative
950 // value for `local_row` and can ask for the length of the row.
951 if (local_row >= 0)
952 return graph->NumMyIndices(local_row);
953 else
954 return static_cast<size_type>(-1);
955 }
956
957
958
959 void
961 const ArrayView<const size_type> &columns,
962 const bool indices_are_sorted)
963 {
964 add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
965 }
966
967
968
969 const Epetra_Map &
971 {
972 // TODO A dynamic_cast fails here, this is suspicious.
973 const auto &domain_map =
974 static_cast<const Epetra_Map &>(graph->DomainMap()); // NOLINT
975 return domain_map;
976 }
977
978
979
980 const Epetra_Map &
982 {
983 // TODO A dynamic_cast fails here, this is suspicious.
984 const auto &range_map =
985 static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
986 return range_map;
987 }
988
989
990
993 {
994 const Epetra_MpiComm *mpi_comm =
995 dynamic_cast<const Epetra_MpiComm *>(&graph->RangeMap().Comm());
996 Assert(mpi_comm != nullptr, ExcInternalError());
997 return mpi_comm->Comm();
998 }
999
1000
1001
1002 void
1007
1008
1009
1010 // As of now, no particularly neat
1011 // output is generated in case of
1012 // multiple processors.
1013 void
1014 SparsityPattern::print(std::ostream &out,
1015 const bool write_extended_trilinos_info) const
1016 {
1017 if (write_extended_trilinos_info)
1018 out << *graph;
1019 else
1020 {
1021 int *indices;
1022 int num_entries;
1023
1024 for (int i = 0; i < graph->NumMyRows(); ++i)
1025 {
1026 graph->ExtractMyRowView(i, num_entries, indices);
1027 for (int j = 0; j < num_entries; ++j)
1028 out << "(" << TrilinosWrappers::global_index(graph->RowMap(), i)
1029 << ","
1030 << TrilinosWrappers::global_index(graph->ColMap(), indices[j])
1031 << ") " << std::endl;
1032 }
1033 }
1034
1035 AssertThrow(out.fail() == false, ExcIO());
1036 }
1037
1038
1039
1040 void
1041 SparsityPattern::print_gnuplot(std::ostream &out) const
1042 {
1043 Assert(graph->Filled() == true, ExcInternalError());
1044 for (::types::global_dof_index row = 0; row < local_size(); ++row)
1045 {
1046 int *indices;
1047 int num_entries;
1048 graph->ExtractMyRowView(row, num_entries, indices);
1049
1050 Assert(num_entries >= 0, ExcInternalError());
1051 // avoid sign comparison warning
1052 const ::types::global_dof_index num_entries_ = num_entries;
1053 for (::types::global_dof_index j = 0; j < num_entries_; ++j)
1054 // while matrix entries are usually
1055 // written (i,j), with i vertical and
1056 // j horizontal, gnuplot output is
1057 // x-y, that is we have to exchange
1058 // the order of output
1059 out << static_cast<int>(
1060 TrilinosWrappers::global_index(graph->ColMap(), indices[j]))
1061 << " "
1062 << -static_cast<int>(
1063 TrilinosWrappers::global_index(graph->RowMap(), row))
1064 << std::endl;
1065 }
1066
1067 AssertThrow(out.fail() == false, ExcIO());
1068 }
1069
1070 // TODO: Implement!
1071 std::size_t
1073 {
1075 return 0;
1076 }
1077
1078
1079# ifndef DOXYGEN
1080 // explicit instantiations
1081 //
1082 template void
1083 SparsityPattern::copy_from(const ::SparsityPattern &);
1084 template void
1085 SparsityPattern::copy_from(const ::DynamicSparsityPattern &);
1086
1087 template void
1089 const ::SparsityPattern &,
1090 const MPI_Comm,
1091 bool);
1092 template void
1094 const ::DynamicSparsityPattern &,
1095 const MPI_Comm,
1096 bool);
1097
1098
1099 template void
1101 const IndexSet &,
1102 const ::SparsityPattern &,
1103 const MPI_Comm,
1104 bool);
1105 template void
1107 const IndexSet &,
1108 const ::DynamicSparsityPattern &,
1109 const MPI_Comm,
1110 bool);
1111# endif
1112
1113} // namespace TrilinosWrappers
1114
1116
1117#endif // DEAL_II_WITH_TRILINOS
iterator begin() const
Definition array_view.h:702
iterator end() const
Definition array_view.h:711
size_type size() const
Definition index_set.h:1776
size_type n_elements() const
Definition index_set.h:1934
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
void subtract_set(const IndexSet &other)
Definition index_set.cc:473
virtual void resize(const size_type rows, const size_type cols)
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
std::pair< size_type, size_type > local_range() const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
void copy_from(const SparsityPattern &input_sparsity_pattern)
std::unique_ptr< Epetra_Map > column_space_map
const_iterator begin() const
bool in_local_range(const size_type index) 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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
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)
#define DEAL_II_NOT_IMPLEMENTED()
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1204
void reinit_sp(const Teuchos::RCP< TpetraTypes::MapType< MemorySpace > > &row_map, const Teuchos::RCP< TpetraTypes::MapType< MemorySpace > > &col_map, const size_type< MemorySpace > n_entries_per_row, Teuchos::RCP< TpetraTypes::MapType< MemorySpace > > &column_space_map, Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > &graph, Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > &nonlocal_graph)
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:92
const Epetra_Comm & comm_self()
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition types.h:32