Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2836-g362046b457 2025-03-13 15:40:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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 if constexpr (running_in_debug_mode())
561 {
562 {
563 IndexSet tmp = writable_rows & row_parallel_partitioning;
564 Assert(tmp == row_parallel_partitioning,
566 "The set of writable rows passed to this method does not "
567 "contain the locally owned rows, which is not allowed."));
568 }
569 }
570 nonlocal_partitioner.subtract_set(row_parallel_partitioning);
571 if (Utilities::MPI::n_mpi_processes(communicator) > 1)
572 {
573 Epetra_Map nonlocal_map =
574 nonlocal_partitioner.make_trilinos_map(communicator, true);
576 std::make_unique<Epetra_CrsGraph>(Copy, nonlocal_map, 0);
577 }
578 else
579 Assert(nonlocal_partitioner.n_elements() == 0, ExcInternalError());
580 }
581
582
583
584 template <typename SparsityPatternType>
585 void
587 const IndexSet &row_parallel_partitioning,
588 const IndexSet &col_parallel_partitioning,
589 const SparsityPatternType &nontrilinos_sparsity_pattern,
590 const MPI_Comm communicator,
591 const bool exchange_data)
592 {
593 SparsityPatternBase::resize(row_parallel_partitioning.size(),
594 col_parallel_partitioning.size());
595 Epetra_Map row_map =
596 row_parallel_partitioning.make_trilinos_map(communicator, false);
597 Epetra_Map col_map =
598 col_parallel_partitioning.make_trilinos_map(communicator, false);
599 reinit_sp(row_map,
600 col_map,
601 nontrilinos_sparsity_pattern,
602 exchange_data,
604 graph,
606 }
607
608
609
610 template <typename SparsityPatternType>
611 void
613 const IndexSet &parallel_partitioning,
614 const SparsityPatternType &nontrilinos_sparsity_pattern,
615 const MPI_Comm communicator,
616 const bool exchange_data)
617 {
618 AssertDimension(nontrilinos_sparsity_pattern.n_rows(),
619 parallel_partitioning.size());
620 AssertDimension(nontrilinos_sparsity_pattern.n_cols(),
621 parallel_partitioning.size());
622 SparsityPatternBase::resize(parallel_partitioning.size(),
623 parallel_partitioning.size());
624 Epetra_Map map =
625 parallel_partitioning.make_trilinos_map(communicator, false);
626 reinit_sp(map,
627 map,
628 nontrilinos_sparsity_pattern,
629 exchange_data,
631 graph,
633 }
634
635
636
639 {
641 return *this;
642 }
643
644
645
646 void
648 {
650 column_space_map = std::make_unique<Epetra_Map>(*sp.column_space_map);
651 graph = std::make_unique<Epetra_FECrsGraph>(*sp.graph);
652
653 if (sp.nonlocal_graph.get() != nullptr)
654 nonlocal_graph = std::make_unique<Epetra_CrsGraph>(*sp.nonlocal_graph);
655 else
656 nonlocal_graph.reset();
657 }
658
659
660
661 template <typename SparsityPatternType>
662 void
663 SparsityPattern::copy_from(const SparsityPatternType &sp)
664 {
665 SparsityPatternBase::resize(sp.n_rows(), sp.n_cols());
666 const Epetra_Map rows(TrilinosWrappers::types::int_type(sp.n_rows()),
667 0,
669 const Epetra_Map columns(TrilinosWrappers::types::int_type(sp.n_cols()),
670 0,
672
673 reinit_sp(
674 rows, columns, sp, false, column_space_map, graph, nonlocal_graph);
675 }
676
677
678
679 void
681 {
683 // When we clear the matrix, reset
684 // the pointer and generate an
685 // empty sparsity pattern.
687 std::make_unique<Epetra_Map>(TrilinosWrappers::types::int_type(0),
690 graph = std::make_unique<Epetra_FECrsGraph>(View,
693 0);
694 graph->FillComplete();
695
696 nonlocal_graph.reset();
697 }
698
699
700
701 void
703 {
704 int ierr;
706 if (nonlocal_graph.get() != nullptr)
707 {
708 if (nonlocal_graph->IndicesAreGlobal() == false &&
709 nonlocal_graph->RowMap().NumMyElements() > 0 &&
711 {
712 // Insert dummy element at (row, column) that corresponds to row 0
713 // in local index counting.
717
718 // in case we have a square sparsity pattern, add the entry on the
719 // diagonal
722 column = row;
723 // if not, take a column index that we have ourselves since we
724 // know for sure it is there (and it will not create spurious
725 // messages to many ranks like putting index 0 on many processors)
726 else if (column_space_map->NumMyElements() > 0)
728 ierr = nonlocal_graph->InsertGlobalIndices(row, 1, &column);
729 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
730 }
731 Assert(nonlocal_graph->RowMap().NumMyElements() == 0 ||
733 nonlocal_graph->IndicesAreGlobal() == true,
735
736 ierr =
737 nonlocal_graph->FillComplete(*column_space_map, graph->RangeMap());
738 AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
739 ierr = nonlocal_graph->OptimizeStorage();
740 AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
741 Epetra_Export exporter(nonlocal_graph->RowMap(), graph->RowMap());
742 ierr = graph->Export(*nonlocal_graph, exporter, Add);
743 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
744 ierr = graph->FillComplete(*column_space_map, graph->RangeMap());
745 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
746 }
747 else
748 {
749 // TODO A dynamic_cast fails here, this is suspicious.
750 const auto &range_map =
751 static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
752 ierr = graph->GlobalAssemble(*column_space_map, range_map, true);
753 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
754 }
755
756 try
757 {
758 ierr = graph->OptimizeStorage();
759 }
760 catch (const int error_code)
761 {
763 false,
765 "The Epetra_CrsGraph::OptimizeStorage() function "
766 "has thrown an error with code " +
767 std::to_string(error_code) +
768 ". You will have to look up the exact meaning of this error "
769 "in the Trilinos source code, but oftentimes, this function "
770 "throwing an error indicates that you are trying to allocate "
771 "more than 2,147,483,647 nonzero entries in the sparsity "
772 "pattern on the local process; this will not work because "
773 "Epetra indexes entries with a simple 'signed int'."));
774 }
775 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
776
777 // Check consistency between the sizes set at the beginning and what
778 // Trilinos stores:
779 using namespace deal_II_exceptions::internals;
780 Assert(compare_for_equality(n_rows(), n_global_rows(*graph)),
782 Assert(compare_for_equality(n_cols(), n_global_cols(*graph)),
784 }
785
786
787
788 bool
790 {
791 return graph->RowMap().LID(
792 static_cast<TrilinosWrappers::types::int_type>(i)) != -1;
793 }
794
795
796
797 bool
799 {
800 if (!row_is_stored_locally(i))
801 {
802 return false;
803 }
804 else
805 {
806 // Extract local indices in
807 // the matrix.
808 int trilinos_i =
809 graph->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
810 trilinos_j =
811 graph->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
812
813 // Check whether the matrix
814 // already is transformed to
815 // local indices.
816 if (graph->Filled() == false)
817 {
818 int nnz_present = graph->NumGlobalIndices(i);
819 int nnz_extracted;
821
822 // Generate the view and make
823 // sure that we have not generated
824 // an error.
825 // TODO: trilinos_i is the local row index -> it is an int but
826 // ExtractGlobalRowView requires trilinos_i to be the global row
827 // index and thus it should be a long long int
828 int ierr = graph->ExtractGlobalRowView(trilinos_i,
829 nnz_extracted,
830 col_indices);
831 (void)ierr;
832 Assert(ierr == 0, ExcTrilinosError(ierr));
833 Assert(nnz_present == nnz_extracted,
834 ExcDimensionMismatch(nnz_present, nnz_extracted));
835
836 // Search the index
837 const std::ptrdiff_t local_col_index =
838 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
839 col_indices;
840
841 if (local_col_index == nnz_present)
842 return false;
843 }
844 else
845 {
846 // Prepare pointers for extraction
847 // of a view of the row.
848 int nnz_present = graph->NumGlobalIndices(i);
849 int nnz_extracted;
850 int *col_indices;
851
852 // Generate the view and make
853 // sure that we have not generated
854 // an error.
855 int ierr =
856 graph->ExtractMyRowView(trilinos_i, nnz_extracted, col_indices);
857 (void)ierr;
858 Assert(ierr == 0, ExcTrilinosError(ierr));
859
860 Assert(nnz_present == nnz_extracted,
861 ExcDimensionMismatch(nnz_present, nnz_extracted));
862
863 // Search the index
864 const std::ptrdiff_t local_col_index =
865 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
866 col_indices;
867
868 if (local_col_index == nnz_present)
869 return false;
870 }
871 }
872
873 return true;
874 }
875
876
877
880 {
881 size_type local_b = 0;
882 for (int i = 0; i < static_cast<int>(local_size()); ++i)
883 {
884 int *indices;
885 int num_entries;
886 graph->ExtractMyRowView(i, num_entries, indices);
887 for (unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
888 ++j)
889 {
890 if (static_cast<size_type>(std::abs(i - indices[j])) > local_b)
891 local_b = std::abs(i - indices[j]);
892 }
893 }
894
896 graph->Comm().MaxAll(reinterpret_cast<TrilinosWrappers::types::int_type *>(
897 &local_b),
898 &global_b,
899 1);
900 return static_cast<size_type>(global_b);
901 }
902
903
904
905 unsigned int
907 {
908 return graph->NumMyRows();
909 }
910
911
912
913 std::pair<SparsityPattern::size_type, SparsityPattern::size_type>
915 {
917 const size_type end = TrilinosWrappers::max_my_gid(graph->RowMap()) + 1;
918
919 return {begin, end};
920 }
921
922
923
924 std::uint64_t
929
930
931
932 unsigned int
934 {
935 return graph->MaxNumIndices();
936 }
937
938
939
942 {
943 Assert(row < n_rows(), ExcInternalError());
944
945 // Get a representation of the where the present row is located on
946 // the current processor
948 graph->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
949
950 // On the processor who owns this row, we'll have a non-negative
951 // value for `local_row` and can ask for the length of the row.
952 if (local_row >= 0)
953 return graph->NumMyIndices(local_row);
954 else
955 return static_cast<size_type>(-1);
956 }
957
958
959
960 void
962 const ArrayView<const size_type> &columns,
963 const bool indices_are_sorted)
964 {
965 add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
966 }
967
968
969
970 const Epetra_Map &
972 {
973 // TODO A dynamic_cast fails here, this is suspicious.
974 const auto &domain_map =
975 static_cast<const Epetra_Map &>(graph->DomainMap()); // NOLINT
976 return domain_map;
977 }
978
979
980
981 const Epetra_Map &
983 {
984 // TODO A dynamic_cast fails here, this is suspicious.
985 const auto &range_map =
986 static_cast<const Epetra_Map &>(graph->RangeMap()); // NOLINT
987 return range_map;
988 }
989
990
991
994 {
995 const Epetra_MpiComm *mpi_comm =
996 dynamic_cast<const Epetra_MpiComm *>(&graph->RangeMap().Comm());
997 Assert(mpi_comm != nullptr, ExcInternalError());
998 return mpi_comm->Comm();
999 }
1000
1001
1002
1003 void
1008
1009
1010
1011 // As of now, no particularly neat
1012 // output is generated in case of
1013 // multiple processors.
1014 void
1015 SparsityPattern::print(std::ostream &out,
1016 const bool write_extended_trilinos_info) const
1017 {
1018 if (write_extended_trilinos_info)
1019 out << *graph;
1020 else
1021 {
1022 int *indices;
1023 int num_entries;
1024
1025 for (int i = 0; i < graph->NumMyRows(); ++i)
1026 {
1027 graph->ExtractMyRowView(i, num_entries, indices);
1028 for (int j = 0; j < num_entries; ++j)
1029 out << "(" << TrilinosWrappers::global_index(graph->RowMap(), i)
1030 << ","
1031 << TrilinosWrappers::global_index(graph->ColMap(), indices[j])
1032 << ") " << std::endl;
1033 }
1034 }
1035
1036 AssertThrow(out.fail() == false, ExcIO());
1037 }
1038
1039
1040
1041 void
1042 SparsityPattern::print_gnuplot(std::ostream &out) const
1043 {
1044 Assert(graph->Filled() == true, ExcInternalError());
1045 for (::types::global_dof_index row = 0; row < local_size(); ++row)
1046 {
1047 int *indices;
1048 int num_entries;
1049 graph->ExtractMyRowView(row, num_entries, indices);
1050
1051 Assert(num_entries >= 0, ExcInternalError());
1052 // avoid sign comparison warning
1053 const ::types::global_dof_index num_entries_ = num_entries;
1054 for (::types::global_dof_index j = 0; j < num_entries_; ++j)
1055 // while matrix entries are usually
1056 // written (i,j), with i vertical and
1057 // j horizontal, gnuplot output is
1058 // x-y, that is we have to exchange
1059 // the order of output
1060 out << static_cast<int>(
1061 TrilinosWrappers::global_index(graph->ColMap(), indices[j]))
1062 << " "
1063 << -static_cast<int>(
1064 TrilinosWrappers::global_index(graph->RowMap(), row))
1065 << std::endl;
1066 }
1067
1068 AssertThrow(out.fail() == false, ExcIO());
1069 }
1070
1071 // TODO: Implement!
1072 std::size_t
1074 {
1076 return 0;
1077 }
1078
1079
1080# ifndef DOXYGEN
1081 // explicit instantiations
1082 //
1083 template void
1084 SparsityPattern::copy_from(const ::SparsityPattern &);
1085 template void
1086 SparsityPattern::copy_from(const ::DynamicSparsityPattern &);
1087
1088 template void
1090 const ::SparsityPattern &,
1091 const MPI_Comm,
1092 bool);
1093 template void
1095 const ::DynamicSparsityPattern &,
1096 const MPI_Comm,
1097 bool);
1098
1099
1100 template void
1102 const IndexSet &,
1103 const ::SparsityPattern &,
1104 const MPI_Comm,
1105 bool);
1106 template void
1108 const IndexSet &,
1109 const ::DynamicSparsityPattern &,
1110 const MPI_Comm,
1111 bool);
1112# endif
1113
1114} // namespace TrilinosWrappers
1115
1117
1118#endif // DEAL_II_WITH_TRILINOS
iterator begin() const
Definition array_view.h:707
iterator end() const
Definition array_view.h:716
size_type size() const
Definition index_set.h:1787
size_type n_elements() const
Definition index_set.h:1945
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:478
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:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_NOT_IMPLEMENTED()
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)
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1215
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:99
const Epetra_Comm & comm_self()
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition types.h:32