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_tpetra_sparsity_pattern.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 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
15#include <deal.II/base/config.h>
16
17#ifdef DEAL_II_TRILINOS_WITH_TPETRA
18
19# include <deal.II/base/mpi.h>
21
26
27# include <limits>
28
30
31namespace LinearAlgebra
32{
33
34 namespace TpetraWrappers
35 {
37 {
38 template <typename MemorySpace>
39 void
41 {
42 // if we are asked to visit the past-the-end line, then simply
43 // release all our caches and go on with life
44 if (static_cast<std::size_t>(this->a_row) == sparsity_pattern->n_rows())
45 {
46 colnum_cache.reset();
47 return;
48 }
49
50 // otherwise first flush Trilinos caches if necessary
51 if (!sparsity_pattern->is_compressed())
52 sparsity_pattern->compress();
53
54 colnum_cache =
55 std::make_shared<std::vector<::types::signed_global_dof_index>>(
56 sparsity_pattern->row_length(this->a_row));
57
58 if (colnum_cache->size() > 0)
59 {
60 // get a representation of the present row
61 std::size_t ncols;
62# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
64 MemorySpace>::nonconst_global_inds_host_view_type
65 column_indices_view(colnum_cache->data(), colnum_cache->size());
66# else
67 Teuchos::ArrayView<::types::signed_global_dof_index>
68 column_indices_view(colnum_cache->data(), colnum_cache->size());
69# endif
70 sparsity_pattern->graph->getGlobalRowCopy(this->a_row,
71 column_indices_view,
72 ncols);
73 AssertThrow(ncols == colnum_cache->size(), ExcInternalError());
74 }
75 }
76 } // namespace SparsityPatternIterators
77
78
79 // The constructor is actually the only point where we have to check whether
80 // we build a serial or a parallel Trilinos matrix. Actually, it does not
81 // even matter how many threads there are, but only if we use an MPI
82 // compiler or a standard compiler. So, even one thread on a configuration
83 // with MPI will still get a parallel interface.
84 template <typename MemorySpace>
98
99
100
101 template <typename MemorySpace>
103 const size_type m,
104 const size_type n,
105 const size_type n_entries_per_row)
106 {
107 reinit(m, n, n_entries_per_row);
108 }
109
110
111
112 template <typename MemorySpace>
114 const size_type m,
115 const size_type n,
116 const std::vector<size_type> &n_entries_per_row)
117 {
118 reinit(m, n, n_entries_per_row);
119 }
120
121
122
123 template <typename MemorySpace>
125 SparsityPattern<MemorySpace> &&other) noexcept
126 : SparsityPatternBase(std::move(other))
127 , column_space_map(std::move(other.column_space_map))
128 , graph(std::move(other.graph))
129 , nonlocal_graph(std::move(other.nonlocal_graph))
130 {}
131
132
133
134 // Copy function only works if the sparsity pattern is empty.
135 template <typename MemorySpace>
137 const SparsityPattern<MemorySpace> &input_sparsity)
138 : SparsityPatternBase(input_sparsity)
139 , column_space_map(Utilities::Trilinos::internal::make_rcp<
140 TpetraTypes::MapType<MemorySpace>>(
141 0,
142 0,
143 Utilities::Trilinos::tpetra_comm_self()))
144 , graph(Utilities::Trilinos::internal::make_rcp<
145 TpetraTypes::GraphType<MemorySpace>>(column_space_map,
146 column_space_map,
147 0))
148 {
149 (void)input_sparsity;
150 Assert(input_sparsity.n_rows() == 0,
152 "Copy constructor only works for empty sparsity patterns."));
153 }
154
155
156
157 template <typename MemorySpace>
159 const IndexSet &parallel_partitioning,
160 const MPI_Comm communicator,
161 const size_type n_entries_per_row)
162 {
163 reinit(parallel_partitioning,
164 parallel_partitioning,
165 communicator,
166 n_entries_per_row);
167 }
168
169
171 template <typename MemorySpace>
173 const IndexSet &parallel_partitioning,
174 const MPI_Comm communicator,
175 const std::vector<size_type> &n_entries_per_row)
176 {
177 reinit(parallel_partitioning,
178 parallel_partitioning,
179 communicator,
180 n_entries_per_row);
181 }
182
183
184
185 template <typename MemorySpace>
187 const IndexSet &row_parallel_partitioning,
188 const IndexSet &col_parallel_partitioning,
189 const MPI_Comm communicator,
190 const size_type n_entries_per_row)
191 {
192 reinit(row_parallel_partitioning,
193 col_parallel_partitioning,
194 communicator,
195 n_entries_per_row);
196 }
197
198
199
200 template <typename MemorySpace>
202 const IndexSet &row_parallel_partitioning,
203 const IndexSet &col_parallel_partitioning,
204 const MPI_Comm communicator,
205 const std::vector<size_type> &n_entries_per_row)
206 {
207 reinit(row_parallel_partitioning,
208 col_parallel_partitioning,
209 communicator,
210 n_entries_per_row);
211 }
212
213
214
215 template <typename MemorySpace>
217 const IndexSet &row_parallel_partitioning,
218 const IndexSet &col_parallel_partitioning,
219 const IndexSet &writable_rows,
220 const MPI_Comm communicator,
221 const size_type n_max_entries_per_row)
222 {
223 reinit(row_parallel_partitioning,
224 col_parallel_partitioning,
225 writable_rows,
226 communicator,
227 n_max_entries_per_row);
228 }
229
230
231
232 template <typename MemorySpace>
233 void
235 const size_type n,
236 const size_type n_entries_per_row)
237 {
238 reinit(complete_index_set(m),
240 MPI_COMM_SELF,
241 n_entries_per_row);
242 }
243
244
245
246 template <typename MemorySpace>
247 void
249 const size_type m,
250 const size_type n,
251 const std::vector<size_type> &n_entries_per_row)
252 {
253 reinit(complete_index_set(m),
255 MPI_COMM_SELF,
256 n_entries_per_row);
257 }
258
259
260
261 namespace SparsityPatternImpl
262 {
263 template <typename MemorySpace>
265
266 template <typename MemorySpace>
267 void
269 const Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> &row_map,
270 const Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> &col_map,
271 const size_type<MemorySpace> n_entries_per_row,
272 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> &column_space_map,
273 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> &graph,
274 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> &nonlocal_graph)
275 {
276 Assert(row_map->isOneToOne(),
277 ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
278 "the maps of different processors."));
279 Assert(col_map->isOneToOne(),
280 ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
281 "the maps of different processors."));
282
283 nonlocal_graph.reset();
284 graph.reset();
285 column_space_map = col_map;
286
287 // for more than one processor, need to specify only row map first and
288 // let the matrix entries decide about the column map (which says which
289 // columns are present in the matrix, not to be confused with the
290 // col_map that tells how the domain dofs of the matrix will be
291 // distributed). for only one processor, we can directly assign the
292 // columns as well. If we use a recent Trilinos version, we can also
293 // require building a non-local graph which gives us thread-safe
294 // initialization.
297 row_map,
298 n_entries_per_row);
299 }
300
301
302
303 template <typename MemorySpace>
304 void
306 const Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> &row_map,
307 const Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> &col_map,
308 const std::vector<size_type<MemorySpace>> &n_entries_per_row,
309 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> &column_space_map,
310 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> &graph,
311 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> &nonlocal_graph)
312 {
313 Assert(row_map->isOneToOne(),
314 ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
315 "the maps of different processors."));
316 Assert(col_map->isOneToOne(),
317 ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
318 "the maps of different processors."));
319
320 // release memory before reallocation
321 nonlocal_graph.reset();
322 graph.reset();
323 AssertDimension(n_entries_per_row.size(),
324 row_map->getGlobalNumElements());
325
326 column_space_map = col_map;
327
328 // Translate the vector of row lengths into one that only stores
329 // those entries that related to the locally stored rows of the matrix:
330 Kokkos::DualView<size_t *, typename MemorySpace::kokkos_space>
331 local_entries_per_row("local_entries_per_row",
332 row_map->getMaxGlobalIndex() -
333 row_map->getMinGlobalIndex());
335 auto local_entries_per_row_host =
336 local_entries_per_row
337 .template view<Kokkos::DefaultHostExecutionSpace>();
338
339 std::uint64_t total_size = 0;
340 for (unsigned int i = 0; i < local_entries_per_row.extent(0); ++i)
341 {
342 local_entries_per_row_host(i) =
343 n_entries_per_row[row_map->getMinGlobalIndex() + i];
344 total_size += local_entries_per_row_host[i];
345 }
346 local_entries_per_row
347 .template modify<Kokkos::DefaultHostExecutionSpace>();
348 local_entries_per_row
349 .template sync<typename MemorySpace::kokkos_space>();
350
352 total_size < static_cast<std::uint64_t>(
353 std::numeric_limits<
356 "You are requesting to store more elements than global ordinal type allows."));
357
360 col_map,
361 local_entries_per_row);
362 }
363
364
365
366 template <typename SparsityPatternType, typename MemorySpace>
367 void
369 const Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> &row_map,
370 const Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> &col_map,
371 const SparsityPatternType &sp,
372 [[maybe_unused]] const bool exchange_data,
373 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> &column_space_map,
374 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> &graph,
375 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> &nonlocal_graph)
376 {
377 nonlocal_graph.reset();
378 graph.reset();
379
380 AssertDimension(sp.n_rows(), row_map->getGlobalNumElements());
381 AssertDimension(sp.n_cols(), col_map->getGlobalNumElements());
382
385
386 Assert(row_map->isContiguous() == true,
388 "This function only works if the row map is contiguous."));
389
390 const size_type<MemorySpace> first_row = row_map->getMinGlobalIndex(),
391 last_row =
392 row_map->getMaxGlobalIndex() + 1;
393 Teuchos::Array<size_t> n_entries_per_row(last_row - first_row);
394
395 for (size_type<MemorySpace> row = first_row; row < last_row; ++row)
396 n_entries_per_row[row - first_row] = sp.row_length(row);
397
399 std::accumulate(n_entries_per_row.begin(),
400 n_entries_per_row.end(),
401 std::uint64_t(0)) <
402 static_cast<std::uint64_t>(std::numeric_limits<int>::max()),
404 "The TrilinosWrappers use Tpetra internally, and "
405 "Trilinos/Tpetra was compiled with 'local ordinate = int'. "
406 "Therefore, 'signed int' is used to represent local indices, "
407 "and only 2,147,483,647 nonzero matrix entries can be stored "
408 "on a single process, but you are requesting more than "
409 "that. Either use more MPI processes or recompile Trilinos "
410 "with 'local ordinate = long long' "));
411
412# if DEAL_II_TRILINOS_VERSION_GTE(12, 16, 0)
413 if (row_map->getComm()->getSize() > 1)
415 TpetraTypes::GraphType<MemorySpace>>(row_map, n_entries_per_row);
416 else
419 col_map,
420 n_entries_per_row);
421# else
422 if (row_map->getComm()->getSize() > 1)
425 row_map, Teuchos::arcpFromArray(n_entries_per_row));
426 else
429 row_map, col_map, Teuchos::arcpFromArray(n_entries_per_row));
431# endif
432
433 AssertDimension(sp.n_rows(), graph->getGlobalNumRows());
434 AssertDimension(sp.n_cols(), graph->getGlobalNumEntries());
435
436 std::vector<TrilinosWrappers::types::int_type> row_indices;
437
438 for (size_type<MemorySpace> row = first_row; row < last_row; ++row)
439 {
440 const TrilinosWrappers::types::int_type row_length =
441 sp.row_length(row);
442 if (row_length == 0)
443 continue;
444
445 row_indices.resize(row_length, -1);
446 {
447 typename SparsityPatternType::iterator p = sp.begin(row);
448 // avoid incrementing p over the end of the current row because
449 // it is slow for DynamicSparsityPattern in parallel
450 for (int col = 0; col < row_length;)
451 {
452 row_indices[col++] = p->column();
453 if (col < row_length)
454 ++p;
455 }
456 }
457 graph->insertGlobalIndices(row, row_length, row_indices.data());
458 }
459
460 graph->globalAssemble();
461 }
462 } // namespace SparsityPatternImpl
464
465 template <typename MemorySpace>
466 void
467 SparsityPattern<MemorySpace>::reinit(const IndexSet &parallel_partitioning,
468 const MPI_Comm communicator,
469 const size_type n_entries_per_row)
470 {
471 SparsityPatternBase::resize(parallel_partitioning.size(),
472 parallel_partitioning.size());
473 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map =
474 parallel_partitioning
475 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
476 communicator, false);
477 SparsityPatternImpl::reinit_sp<MemorySpace>(
478 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
479 }
480
482
483 template <typename MemorySpace>
484 void
486 const IndexSet &parallel_partitioning,
487 const MPI_Comm communicator,
488 const std::vector<size_type> &n_entries_per_row)
489 {
490 SparsityPatternBase::resize(parallel_partitioning.size(),
491 parallel_partitioning.size());
492 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map =
493 parallel_partitioning
494 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
495 communicator, false);
496 SparsityPatternImpl::reinit_sp<MemorySpace>(
497 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
498 }
499
500
501
502 template <typename MemorySpace>
503 void
505 const IndexSet &row_parallel_partitioning,
506 const IndexSet &col_parallel_partitioning,
507 const MPI_Comm communicator,
508 const size_type n_entries_per_row)
509 {
510 SparsityPatternBase::resize(row_parallel_partitioning.size(),
511 col_parallel_partitioning.size());
512 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
513 row_parallel_partitioning
514 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
515 communicator, false);
516 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
517 col_parallel_partitioning
518 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
519 communicator, false);
520 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
521 col_map,
522 n_entries_per_row,
523 column_space_map,
524 graph,
525 nonlocal_graph);
526 }
527
529
530 template <typename MemorySpace>
531 void
533 const IndexSet &row_parallel_partitioning,
534 const IndexSet &col_parallel_partitioning,
535 const MPI_Comm communicator,
536 const std::vector<size_type> &n_entries_per_row)
537 {
538 SparsityPatternBase::resize(row_parallel_partitioning.size(),
539 col_parallel_partitioning.size());
540 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
541 row_parallel_partitioning
542 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
543 communicator, false);
544 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
545 col_parallel_partitioning
546 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
547 communicator, false);
548 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
549 col_map,
550 n_entries_per_row,
551 column_space_map,
552 graph,
553 nonlocal_graph);
554 }
555
556
557
558 template <typename MemorySpace>
559 void
561 const IndexSet &row_parallel_partitioning,
562 const IndexSet &col_parallel_partitioning,
563 const IndexSet &writable_rows,
564 const MPI_Comm communicator,
565 const size_type n_entries_per_row)
566 {
567 SparsityPatternBase::resize(row_parallel_partitioning.size(),
568 col_parallel_partitioning.size());
569 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
570 row_parallel_partitioning
571 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
572 communicator, false);
573 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
574 col_parallel_partitioning
575 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
576 communicator, false);
577 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
578 col_map,
579 n_entries_per_row,
580 column_space_map,
581 graph,
582 nonlocal_graph);
583
584 IndexSet nonlocal_partitioner = writable_rows;
585 AssertDimension(nonlocal_partitioner.size(),
586 row_parallel_partitioning.size());
587 if constexpr (running_in_debug_mode())
588 {
589 {
590 IndexSet tmp = writable_rows & row_parallel_partitioning;
591 Assert(tmp == row_parallel_partitioning,
593 "The set of writable rows passed to this method does not "
594 "contain the locally owned rows, which is not allowed."));
595 }
596 }
597 nonlocal_partitioner.subtract_set(row_parallel_partitioning);
598 if (Utilities::MPI::n_mpi_processes(communicator) > 1)
599 {
600 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> nonlocal_map =
601 nonlocal_partitioner
602 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
603 communicator, true);
605 TpetraTypes::GraphType<MemorySpace>>(nonlocal_map, col_map, 0);
606 }
607 else
608 Assert(nonlocal_partitioner.n_elements() == 0, ExcInternalError());
609 }
610
611
612
613 template <typename MemorySpace>
614 template <typename SparsityPatternType>
615 void
617 const IndexSet &row_parallel_partitioning,
618 const IndexSet &col_parallel_partitioning,
619 const SparsityPatternType &nontrilinos_sparsity_pattern,
620 const MPI_Comm communicator,
621 const bool exchange_data)
622 {
623 SparsityPatternBase::resize(row_parallel_partitioning.size(),
624 col_parallel_partitioning.size());
625 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
626 row_parallel_partitioning
627 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
628 communicator, false);
629 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
630 col_parallel_partitioning
631 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
632 communicator, false);
633 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
634 row_map,
635 col_map,
636 nontrilinos_sparsity_pattern,
637 exchange_data,
638 column_space_map,
639 graph,
640 nonlocal_graph);
641 }
642
643
645 template <typename MemorySpace>
646 template <typename SparsityPatternType>
647 void
649 const IndexSet &parallel_partitioning,
650 const SparsityPatternType &nontrilinos_sparsity_pattern,
651 const MPI_Comm communicator,
652 const bool exchange_data)
653 {
654 AssertDimension(nontrilinos_sparsity_pattern.n_rows(),
655 parallel_partitioning.size());
656 AssertDimension(nontrilinos_sparsity_pattern.n_cols(),
657 parallel_partitioning.size());
658 SparsityPatternBase::resize(parallel_partitioning.size(),
659 parallel_partitioning.size());
660 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map =
661 parallel_partitioning
662 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
663 communicator, false);
664 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
665 map,
666 map,
667 nontrilinos_sparsity_pattern,
668 exchange_data,
669 column_space_map,
670 graph,
671 nonlocal_graph);
672 }
673
674
675
676 template <typename MemorySpace>
680 {
682 return *this;
683 }
684
685
686
687 template <typename MemorySpace>
688 void
691 {
697
698 if (sp.nonlocal_graph.get() != nullptr)
701 else
702 nonlocal_graph.reset();
703 }
705
706
707 template <typename MemorySpace>
708 template <typename SparsityPatternType>
709 void
710 SparsityPattern<MemorySpace>::copy_from(const SparsityPatternType &sp)
711 {
712 SparsityPatternBase::resize(sp.n_rows(), sp.n_cols());
713 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> rows =
717 0,
719 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> columns =
723 0,
725
726 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
727 rows, columns, sp, false, column_space_map, graph, nonlocal_graph);
729
730
731
732 template <typename MemorySpace>
733 void
735 {
737 // When we clear the matrix, reset
738 // the pointer and generate an
739 // empty sparsity pattern.
746 TpetraTypes::GraphType<MemorySpace>>(column_space_map,
747 column_space_map,
748 0);
749 graph->fillComplete();
750
751 nonlocal_graph.reset();
752 }
753
754
755
756 template <typename MemorySpace>
757 void
759 {
760 Assert(column_space_map.get(), ExcInternalError());
761 if (nonlocal_graph.get() != nullptr)
762 {
763# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
764 if (nonlocal_graph->getRowMap()->getLocalNumElements() > 0 &&
765 column_space_map->getGlobalNumElements() > 0)
766# else
767 if (nonlocal_graph->getRowMap()->getNodeNumElements() > 0 &&
768 column_space_map->getGlobalNumElements() > 0)
769# endif
770 {
771 // Insert dummy element at (row, column) that corresponds to row 0
772 // in local index counting.
774 nonlocal_graph->getRowMap()->getGlobalElement(0);
776
777 // in case we have a square sparsity pattern, add the entry on the
778 // diagonal
779 if (column_space_map->getGlobalNumElements() ==
780 graph->getRangeMap()->getGlobalNumElements())
781 column = row;
782 // if not, take a column index that we have ourselves since we
783 // know for sure it is there (and it will not create spurious
784 // messages to many ranks like putting index 0 on many
785 // processors)
786# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
787 else if (column_space_map->getLocalNumElements() > 0)
788# else
789 else if (column_space_map->getNodeNumElements() > 0)
790# endif
791 column = column_space_map->getGlobalElement(0);
792 nonlocal_graph->insertGlobalIndices(row, 1, &column);
793 }
794# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
795 Assert(nonlocal_graph->getRowMap()->getLocalNumElements() == 0 ||
796 column_space_map->getGlobalNumElements() == 0,
798# else
799 Assert(nonlocal_graph->getRowMap()->getNodeNumElements() == 0 ||
800 column_space_map->getGlobalNumElements() == 0,
802# endif
803
804 nonlocal_graph->fillComplete(column_space_map, graph->getRowMap());
805 }
806 graph->fillComplete(column_space_map, graph->getRowMap());
807
808 // Check consistency between the sizes set at the beginning and what
809 // Trilinos stores:
810 using namespace deal_II_exceptions::internals;
811 AssertDimension(n_rows(), graph->getGlobalNumRows());
812 AssertDimension(n_cols(), graph->getGlobalNumCols());
813 }
814
815
816
817 template <typename MemorySpace>
818 bool
820 {
821 return graph->getRowMap()->getLocalElement(i) !=
822 Teuchos::OrdinalTraits<int>::invalid();
823 }
824
825
827 template <typename MemorySpace>
828 bool
830 const size_type j) const
831 {
832# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
833 if (!row_is_stored_locally(i))
834 return false;
835
836 // Extract local indices in the matrix.
837 const auto trilinos_i = graph->getRowMap()->getLocalElement(i);
838 const auto trilinos_j = graph->getColMap()->getLocalElement(j);
839
841 col_indices;
842
843 // Generate the view.
844 graph->getLocalRowView(trilinos_i, col_indices);
845
846 // Search the index
847 const size_type local_col_index =
848 std::find(col_indices.data(),
849 col_indices.data() + col_indices.size(),
850 trilinos_j) -
851 col_indices.data();
852
853 return static_cast<std::size_t>(local_col_index) != col_indices.size();
854# else
855 (void)i;
856 (void)j;
858 return false;
859# endif
860 }
861
862
863
864 template <typename MemorySpace>
867 {
868 size_type local_b = 0;
869 for (int i = 0; i < static_cast<int>(local_size()); ++i)
870 {
871# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
872 typename TpetraTypes::GraphType<
873 MemorySpace>::local_inds_host_view_type indices;
874# else
875 Teuchos::ArrayView<const int> indices;
876# endif
877
878 graph->getLocalRowView(i, indices);
879 const auto num_entries = indices.size();
880 for (unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
881 ++j)
882 {
883 if (static_cast<size_type>(std::abs(i - indices[j])) > local_b)
884 local_b = std::abs(i - indices[j]);
885 }
886 }
887
889 Utilities::MPI::max(local_b,
891 graph->getComm()));
892 return static_cast<size_type>(global_b);
893 }
894
895
896
897 template <typename MemorySpace>
898 unsigned int
900 {
901# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
902 return graph->getLocalNumRows();
903# else
904 return graph->getNodeNumRows();
905# endif
906 }
907
908
909
910 template <typename MemorySpace>
911 std::pair<typename SparsityPattern<MemorySpace>::size_type,
914 {
915 const size_type begin = graph->getRowMap()->getMinGlobalIndex();
916 const size_type end = graph->getRowMap()->getMaxGlobalIndex() + 1;
917
918 return {begin, end};
919 }
920
921
922
923 template <typename MemorySpace>
924 std::uint64_t
926 {
927 return graph->getGlobalNumEntries();
928 }
929
930
931
932 template <typename MemorySpace>
933 unsigned int
935 {
936# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
937 return graph->getLocalMaxNumRowEntries();
938# else
939 return graph->getNodeMaxNumRowEntries();
940# endif
941 }
942
943
944
945 template <typename MemorySpace>
948 {
949 Assert(row < (size_type)n_rows(), ExcInternalError());
950
951 // Get a representation of the where the present row is located on
952 // the current processor
954 graph->getRowMap()->getLocalElement(row);
955
956 // On the processor who owns this row, we'll have a non-negative
957 // value for `local_row` and can ask for the length of the row.
958 if (local_row >= 0)
959 return graph->getNumEntriesInLocalRow(local_row);
960 else
961 return static_cast<size_type>(-1);
962 }
963
964
965
966 template <typename MemorySpace>
967 void
969 const ::types::global_dof_index &row,
971 const bool indices_are_sorted)
972 {
973 add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
974 }
975
976
977
978 template <typename MemorySpace>
979 Teuchos::RCP<const typename TpetraTypes::MapType<MemorySpace>>
981 {
982 return graph->getDomainMap();
983 }
984
985
986
987 template <typename MemorySpace>
988 Teuchos::RCP<const typename TpetraTypes::MapType<MemorySpace>>
990 {
991 return graph->getRangeMap();
992 }
993
994
995
996 template <typename MemorySpace>
999 {
1001 graph->getRangeMap()->getComm());
1002 }
1003
1004
1005
1006 template <typename MemorySpace>
1007 Teuchos::RCP<const Teuchos::Comm<int>>
1009 {
1010 return graph->getRangeMap()->getComm();
1011 }
1012
1013
1014
1015 // As of now, no particularly neat
1016 // output is generated in case of
1017 // multiple processors.
1018 template <typename MemorySpace>
1019 void
1021 std::ostream &out,
1022 const bool write_extended_trilinos_info) const
1023 {
1024 if (write_extended_trilinos_info)
1025 out << *graph;
1026 else
1027 {
1028# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
1029 for (unsigned int i = 0; i < graph->getLocalNumRows(); ++i)
1030# else
1031 for (unsigned int i = 0; i < graph->getNodeNumRows(); ++i)
1032# endif
1033 {
1034# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1035 typename TpetraTypes::GraphType<
1036 MemorySpace>::local_inds_host_view_type indices;
1037# else
1038 Teuchos::ArrayView<const int> indices;
1039# endif
1040 graph->getLocalRowView(i, indices);
1041 int num_entries = indices.size();
1042 for (int j = 0; j < num_entries; ++j)
1043 out << "(" << graph->getRowMap()->getGlobalElement(i) << ","
1044 << graph->getColMap()->getGlobalElement(indices[j]) << ") "
1045 << std::endl;
1046 }
1047 }
1048
1049 AssertThrow(out.fail() == false, ExcIO());
1050 }
1051
1052
1053
1054 template <typename MemorySpace>
1055 void
1057 {
1058 Assert(graph->isFillComplete() == true, ExcInternalError());
1059
1060 for (unsigned int row = 0; row < local_size(); ++row)
1061 {
1062# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1063 typename TpetraTypes::GraphType<
1064 MemorySpace>::local_inds_host_view_type indices;
1065# else
1066 Teuchos::ArrayView<const int> indices;
1067# endif
1068 graph->getLocalRowView(row, indices);
1069 int num_entries = indices.size();
1070
1071 Assert(num_entries >= 0, ExcInternalError());
1072 // avoid sign comparison warning
1073 const ::types::signed_global_dof_index num_entries_ =
1074 num_entries;
1075 for (::types::signed_global_dof_index j = 0; j < num_entries_;
1076 ++j)
1077 // while matrix entries are usually
1078 // written (i,j), with i vertical and
1079 // j horizontal, gnuplot output is
1080 // x-y, that is we have to exchange
1081 // the order of output
1082 out << static_cast<int>(
1083 graph->getColMap()->getGlobalElement(indices[j]))
1084 << " "
1085 << -static_cast<int>(graph->getRowMap()->getGlobalElement(row))
1086 << std::endl;
1087 }
1088
1089 AssertThrow(out.fail() == false, ExcIO());
1090 }
1091
1092 // TODO: Implement!
1093 template <typename MemorySpace>
1094 std::size_t
1100
1101
1102# ifndef DOXYGEN
1103 // explicit instantiations
1105
1106 template void
1108 const ::SparsityPattern &);
1109 template void
1111 const ::DynamicSparsityPattern &);
1112
1113 template void
1115 const IndexSet &,
1116 const ::SparsityPattern &,
1117 const MPI_Comm,
1118 bool);
1119 template void
1121 const IndexSet &,
1122 const ::DynamicSparsityPattern &,
1123 const MPI_Comm,
1124 bool);
1125
1126
1127 template void
1129 const IndexSet &,
1130 const IndexSet &,
1131 const ::SparsityPattern &,
1132 const MPI_Comm,
1133 bool);
1134 template void
1136 const IndexSet &,
1137 const IndexSet &,
1138 const ::DynamicSparsityPattern &,
1139 const MPI_Comm,
1140 bool);
1141
1142
1144
1145 template void
1147 const ::SparsityPattern &);
1148 template void
1150 const ::DynamicSparsityPattern &);
1151
1152 template void
1154 const IndexSet &,
1155 const ::SparsityPattern &,
1156 const MPI_Comm,
1157 bool);
1158 template void
1160 const IndexSet &,
1161 const ::DynamicSparsityPattern &,
1162 const MPI_Comm,
1163 bool);
1164
1165
1166 template void
1168 const IndexSet &,
1169 const IndexSet &,
1170 const ::SparsityPattern &,
1171 const MPI_Comm,
1172 bool);
1173 template void
1175 const IndexSet &,
1176 const IndexSet &,
1177 const ::DynamicSparsityPattern &,
1178 const MPI_Comm,
1179 bool);
1180
1181# endif
1182
1183 } // namespace TpetraWrappers
1184
1185} // namespace LinearAlgebra
1186
1188
1189#endif // DEAL_II_TRILINOS_WITH_TPETRA
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
Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > graph
Teuchos::RCP< const Teuchos::Comm< int > > get_teuchos_mpi_communicator() const
Teuchos::RCP< const TpetraTypes::MapType< MemorySpace > > range_partitioner() const
void copy_from(const SparsityPattern< MemorySpace > &input_sparsity_pattern)
SparsityPattern< MemorySpace > & operator=(const SparsityPattern< MemorySpace > &input_sparsity_pattern)
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > nonlocal_graph
Teuchos::RCP< const TpetraTypes::MapType< MemorySpace > > domain_partitioner() const
bool exists(const size_type i, const size_type j) const
Teuchos::RCP< TpetraTypes::MapType< MemorySpace > > column_space_map
virtual void add_row_entries(const ::types::global_dof_index &row, const ArrayView< const ::types::global_dof_index > &columns, const bool indices_are_sorted=false) override
void reinit(const size_type m, const size_type n, const size_type n_entries_per_row=0)
virtual void resize(const size_type rows, const size_type cols)
size_type n_rows() const
size_type n_cols() const
#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 & 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)
typename SparsityPattern< MemorySpace >::size_type size_type
Tpetra::CrsGraph< LO, GO, NodeType< MemorySpace > > GraphType
Tpetra::Map< LO, GO, NodeType< MemorySpace > > MapType
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:99
T max(const T &t, const MPI_Comm mpi_communicator)
Teuchos::RCP< T > make_rcp(Args &&...args)
MPI_Comm teuchos_comm_to_mpi_comm(const Teuchos::RCP< const Teuchos::Comm< int > > &teuchos_comm)
const Teuchos::RCP< const Teuchos::Comm< int > > & tpetra_comm_self()
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)