Reference documentation for deal.II version GIT relicensing-384-g4769914377 2024-04-13 09:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
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<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)
63 typename Tpetra::CrsGraph<
64 int,
66 NodeType>::nonconst_global_inds_host_view_type
67 column_indices_view(colnum_cache->data(), colnum_cache->size());
68# else
69 Teuchos::ArrayView<::types::signed_global_dof_index>
70 column_indices_view(colnum_cache->data(), colnum_cache->size());
71# endif
72 sparsity_pattern->graph->getGlobalRowCopy(this->a_row,
73 column_indices_view,
74 ncols);
75 AssertThrow(ncols == colnum_cache->size(), ExcInternalError());
76 }
77 }
78 } // namespace SparsityPatternIterators
79
80
81 // The constructor is actually the only point where we have to check whether
82 // we build a serial or a parallel Trilinos matrix. Actually, it does not
83 // even matter how many threads there are, but only if we use an MPI
84 // compiler or a standard compiler. So, even one thread on a configuration
85 // with MPI will still get a parallel interface.
86 template <typename MemorySpace>
88 {
89 column_space_map = Utilities::Trilinos::internal::make_rcp<MapType>(
93 graph =
94 Utilities::Trilinos::internal::make_rcp<GraphType>(column_space_map,
95 column_space_map,
96 0);
97 graph->fillComplete();
98 }
99
100
101
102 template <typename MemorySpace>
104 const size_type m,
105 const size_type n,
106 const size_type n_entries_per_row)
107 {
108 reinit(m, n, n_entries_per_row);
109 }
110
111
112
113 template <typename MemorySpace>
115 const size_type m,
116 const size_type n,
117 const std::vector<size_type> &n_entries_per_row)
118 {
119 reinit(m, n, n_entries_per_row);
120 }
121
122
123
124 template <typename MemorySpace>
126 SparsityPattern<MemorySpace> &&other) noexcept
127 : SparsityPatternBase(std::move(other))
128 , column_space_map(std::move(other.column_space_map))
129 , graph(std::move(other.graph))
130 , nonlocal_graph(std::move(other.nonlocal_graph))
131 {}
132
133
134
135 // Copy function only works if the sparsity pattern is empty.
136 template <typename MemorySpace>
138 const SparsityPattern<MemorySpace> &input_sparsity)
139 : SparsityPatternBase(input_sparsity)
140 , column_space_map(Utilities::Trilinos::internal::make_rcp<MapType>(
141 0,
142 0,
143 Utilities::Trilinos::tpetra_comm_self()))
144 , graph(
145 Utilities::Trilinos::internal::make_rcp<GraphType>(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
170
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 using MapType =
268 Tpetra::Map<int,
271
272 template <typename MemorySpace>
273 using GraphType =
274 Tpetra::CrsGraph<int,
277
278 template <typename MemorySpace>
279 void
280 reinit_sp(const Teuchos::RCP<MapType<MemorySpace>> &row_map,
281 const Teuchos::RCP<MapType<MemorySpace>> &col_map,
282 const size_type<MemorySpace> n_entries_per_row,
283 Teuchos::RCP<MapType<MemorySpace>> &column_space_map,
284 Teuchos::RCP<GraphType<MemorySpace>> &graph,
285 Teuchos::RCP<GraphType<MemorySpace>> &nonlocal_graph)
286 {
287 Assert(row_map->isOneToOne(),
288 ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
289 "the maps of different processors."));
290 Assert(col_map->isOneToOne(),
291 ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
292 "the maps of different processors."));
293
294 nonlocal_graph.reset();
295 graph.reset();
296 column_space_map = col_map;
297
298 // for more than one processor, need to specify only row map first and
299 // let the matrix entries decide about the column map (which says which
300 // columns are present in the matrix, not to be confused with the
301 // col_map that tells how the domain dofs of the matrix will be
302 // distributed). for only one processor, we can directly assign the
303 // columns as well. If we use a recent Trilinos version, we can also
304 // require building a non-local graph which gives us thread-safe
305 // initialization.
306 graph = Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
307 row_map, row_map, n_entries_per_row);
308 }
309
310
311
312 template <typename MemorySpace>
313 void
314 reinit_sp(const Teuchos::RCP<MapType<MemorySpace>> &row_map,
315 const Teuchos::RCP<MapType<MemorySpace>> &col_map,
316 const std::vector<size_type<MemorySpace>> &n_entries_per_row,
317 Teuchos::RCP<MapType<MemorySpace>> &column_space_map,
318 Teuchos::RCP<GraphType<MemorySpace>> &graph,
319 Teuchos::RCP<GraphType<MemorySpace>> &nonlocal_graph)
320 {
321 Assert(row_map->isOneToOne(),
322 ExcMessage("Row map must be 1-to-1, i.e., no overlap between "
323 "the maps of different processors."));
324 Assert(col_map->isOneToOne(),
325 ExcMessage("Column map must be 1-to-1, i.e., no overlap between "
326 "the maps of different processors."));
327
328 // release memory before reallocation
329 nonlocal_graph.reset();
330 graph.reset();
331 AssertDimension(n_entries_per_row.size(),
332 row_map->getGlobalNumElements());
333
334 column_space_map = col_map;
335
336 // Translate the vector of row lengths into one that only stores
337 // those entries that related to the locally stored rows of the matrix:
338 Kokkos::DualView<size_t *, typename MemorySpace::kokkos_space>
339 local_entries_per_row("local_entries_per_row",
340 row_map->getMaxGlobalIndex() -
341 row_map->getMinGlobalIndex());
342
343 auto local_entries_per_row_host =
344 local_entries_per_row
345 .template view<Kokkos::DefaultHostExecutionSpace>();
346
347 std::uint64_t total_size = 0;
348 for (unsigned int i = 0; i < local_entries_per_row.extent(0); ++i)
349 {
350 local_entries_per_row_host(i) =
351 n_entries_per_row[row_map->getMinGlobalIndex() + i];
352 total_size += local_entries_per_row_host[i];
353 }
354 local_entries_per_row
355 .template modify<Kokkos::DefaultHostExecutionSpace>();
356 local_entries_per_row
357 .template sync<typename MemorySpace::kokkos_space>();
358
360 total_size < static_cast<std::uint64_t>(
361 std::numeric_limits<
364 "You are requesting to store more elements than global ordinal type allows."));
365
366 graph = Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
367 row_map, col_map, local_entries_per_row);
368 }
369
370
371
372 template <typename SparsityPatternType, typename MemorySpace>
373 void
374 reinit_sp(const Teuchos::RCP<MapType<MemorySpace>> &row_map,
375 const Teuchos::RCP<MapType<MemorySpace>> &col_map,
376 const SparsityPatternType &sp,
377 [[maybe_unused]] const bool exchange_data,
378 Teuchos::RCP<MapType<MemorySpace>> &column_space_map,
379 Teuchos::RCP<GraphType<MemorySpace>> &graph,
380 Teuchos::RCP<GraphType<MemorySpace>> &nonlocal_graph)
381 {
382 nonlocal_graph.reset();
383 graph.reset();
385 AssertDimension(sp.n_rows(), row_map->getGlobalNumElements());
386 AssertDimension(sp.n_cols(), col_map->getGlobalNumElements());
387
388 column_space_map =
389 Utilities::Trilinos::internal::make_rcp<MapType<MemorySpace>>(
390 *col_map);
391
392 Assert(row_map->isContiguous() == true,
394 "This function only works if the row map is contiguous."));
395
396 const size_type<MemorySpace> first_row = row_map->getMinGlobalIndex(),
397 last_row =
398 row_map->getMaxGlobalIndex() + 1;
399 Teuchos::Array<size_t> n_entries_per_row(last_row - first_row);
400
401 for (size_type<MemorySpace> row = first_row; row < last_row; ++row)
402 n_entries_per_row[row - first_row] = sp.row_length(row);
403
405 std::accumulate(n_entries_per_row.begin(),
406 n_entries_per_row.end(),
407 std::uint64_t(0)) <
408 static_cast<std::uint64_t>(std::numeric_limits<int>::max()),
410 "The TrilinosWrappers use Tpetra internally, and "
411 "Trilinos/Tpetra was compiled with 'local ordinate = int'. "
412 "Therefore, 'signed int' is used to represent local indices, "
413 "and only 2,147,483,647 nonzero matrix entries can be stored "
414 "on a single process, but you are requesting more than "
415 "that. Either use more MPI processes or recompile Trilinos "
416 "with 'local ordinate = long long' "));
417
418# if DEAL_II_TRILINOS_VERSION_GTE(12, 16, 0)
419 if (row_map->getComm()->getSize() > 1)
420 graph =
421 Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
422 row_map, n_entries_per_row);
423 else
424 graph =
425 Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
426 row_map, col_map, n_entries_per_row);
427# else
428 if (row_map->getComm()->getSize() > 1)
429 graph =
431 row_map, Teuchos::arcpFromArray(n_entries_per_row));
432 else
433 graph =
434 Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
435 row_map, col_map, Teuchos::arcpFromArray(n_entries_per_row));
436
437# endif
438
439 AssertDimension(sp.n_rows(), graph->getGlobalNumRows());
440 AssertDimension(sp.n_cols(), graph->getGlobalNumEntries());
441
442 std::vector<TrilinosWrappers::types::int_type> row_indices;
443
444 for (size_type<MemorySpace> row = first_row; row < last_row; ++row)
445 {
446 const TrilinosWrappers::types::int_type row_length =
447 sp.row_length(row);
448 if (row_length == 0)
449 continue;
450
451 row_indices.resize(row_length, -1);
452 {
453 typename SparsityPatternType::iterator p = sp.begin(row);
454 // avoid incrementing p over the end of the current row because
455 // it is slow for DynamicSparsityPattern in parallel
456 for (int col = 0; col < row_length;)
457 {
458 row_indices[col++] = p->column();
459 if (col < row_length)
460 ++p;
461 }
462 }
463 graph->insertGlobalIndices(row, row_length, row_indices.data());
464 }
465
466 graph->globalAssemble();
467 }
468 } // namespace SparsityPatternImpl
469
470
471 template <typename MemorySpace>
472 void
473 SparsityPattern<MemorySpace>::reinit(const IndexSet &parallel_partitioning,
474 const MPI_Comm communicator,
475 const size_type n_entries_per_row)
476 {
477 SparsityPatternBase::resize(parallel_partitioning.size(),
478 parallel_partitioning.size());
479 Teuchos::RCP<MapType> map =
480 parallel_partitioning.make_tpetra_map_rcp(communicator, false);
481 SparsityPatternImpl::reinit_sp<MemorySpace>(
482 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
483 }
484
485
486
487 template <typename MemorySpace>
488 void
490 const IndexSet &parallel_partitioning,
491 const MPI_Comm communicator,
492 const std::vector<size_type> &n_entries_per_row)
493 {
494 SparsityPatternBase::resize(parallel_partitioning.size(),
495 parallel_partitioning.size());
496 Teuchos::RCP<MapType> map =
497 parallel_partitioning.make_tpetra_map_rcp(communicator, false);
498 SparsityPatternImpl::reinit_sp<MemorySpace>(
499 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
500 }
501
502
503
504 template <typename MemorySpace>
505 void
507 const IndexSet &row_parallel_partitioning,
508 const IndexSet &col_parallel_partitioning,
509 const MPI_Comm communicator,
510 const size_type n_entries_per_row)
511 {
512 SparsityPatternBase::resize(row_parallel_partitioning.size(),
513 col_parallel_partitioning.size());
514 Teuchos::RCP<MapType> row_map =
515 row_parallel_partitioning.make_tpetra_map_rcp(communicator, false);
516 Teuchos::RCP<MapType> col_map =
517 col_parallel_partitioning.make_tpetra_map_rcp(communicator, false);
518 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
519 col_map,
520 n_entries_per_row,
521 column_space_map,
522 graph,
523 nonlocal_graph);
524 }
525
526
527
528 template <typename MemorySpace>
529 void
531 const IndexSet &row_parallel_partitioning,
532 const IndexSet &col_parallel_partitioning,
533 const MPI_Comm communicator,
534 const std::vector<size_type> &n_entries_per_row)
535 {
536 SparsityPatternBase::resize(row_parallel_partitioning.size(),
537 col_parallel_partitioning.size());
538 Teuchos::RCP<MapType> row_map =
539 row_parallel_partitioning.make_tpetra_map_rcp(communicator, false);
540 Teuchos::RCP<MapType> col_map =
541 col_parallel_partitioning.make_tpetra_map_rcp(communicator, false);
542 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
543 col_map,
544 n_entries_per_row,
545 column_space_map,
546 graph,
547 nonlocal_graph);
548 }
549
550
551
552 template <typename MemorySpace>
553 void
555 const IndexSet &row_parallel_partitioning,
556 const IndexSet &col_parallel_partitioning,
557 const IndexSet &writable_rows,
558 const MPI_Comm communicator,
559 const size_type n_entries_per_row)
560 {
561 SparsityPatternBase::resize(row_parallel_partitioning.size(),
562 col_parallel_partitioning.size());
563 Teuchos::RCP<MapType> row_map =
564 row_parallel_partitioning.make_tpetra_map_rcp(communicator, false);
565 Teuchos::RCP<MapType> col_map =
566 col_parallel_partitioning.make_tpetra_map_rcp(communicator, false);
567 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
568 col_map,
569 n_entries_per_row,
570 column_space_map,
571 graph,
572 nonlocal_graph);
573
574 IndexSet nonlocal_partitioner = writable_rows;
575 AssertDimension(nonlocal_partitioner.size(),
576 row_parallel_partitioning.size());
577# ifdef DEBUG
579 IndexSet tmp = writable_rows & row_parallel_partitioning;
580 Assert(tmp == row_parallel_partitioning,
582 "The set of writable rows passed to this method does not "
583 "contain the locally owned rows, which is not allowed."));
584 }
585# endif
586 nonlocal_partitioner.subtract_set(row_parallel_partitioning);
587 if (Utilities::MPI::n_mpi_processes(communicator) > 1)
588 {
589 Teuchos::RCP<MapType> nonlocal_map =
590 nonlocal_partitioner.make_tpetra_map_rcp(communicator, true);
591 nonlocal_graph =
592 Utilities::Trilinos::internal::make_rcp<GraphType>(nonlocal_map,
593 col_map,
594 0);
595 }
596 else
597 Assert(nonlocal_partitioner.n_elements() == 0, ExcInternalError());
598 }
599
601
602 template <typename MemorySpace>
603 template <typename SparsityPatternType>
604 void
606 const IndexSet &row_parallel_partitioning,
607 const IndexSet &col_parallel_partitioning,
608 const SparsityPatternType &nontrilinos_sparsity_pattern,
609 const MPI_Comm communicator,
610 const bool exchange_data)
611 {
612 SparsityPatternBase::resize(row_parallel_partitioning.size(),
613 col_parallel_partitioning.size());
614 Teuchos::RCP<MapType> row_map =
615 row_parallel_partitioning.make_tpetra_map_rcp(communicator, false);
616 Teuchos::RCP<MapType> col_map =
617 col_parallel_partitioning.make_tpetra_map_rcp(communicator, false);
618 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
619 row_map,
620 col_map,
621 nontrilinos_sparsity_pattern,
622 exchange_data,
623 column_space_map,
624 graph,
625 nonlocal_graph);
626 }
627
628
629
630 template <typename MemorySpace>
631 template <typename SparsityPatternType>
632 void
634 const IndexSet &parallel_partitioning,
635 const SparsityPatternType &nontrilinos_sparsity_pattern,
636 const MPI_Comm communicator,
637 const bool exchange_data)
638 {
639 AssertDimension(nontrilinos_sparsity_pattern.n_rows(),
640 parallel_partitioning.size());
641 AssertDimension(nontrilinos_sparsity_pattern.n_cols(),
642 parallel_partitioning.size());
643 SparsityPatternBase::resize(parallel_partitioning.size(),
644 parallel_partitioning.size());
645 Teuchos::RCP<MapType> map =
646 parallel_partitioning.make_tpetra_map_rcp(communicator, false);
647 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
648 map,
649 map,
650 nontrilinos_sparsity_pattern,
651 exchange_data,
652 column_space_map,
653 graph,
654 nonlocal_graph);
655 }
656
657
658
659 template <typename MemorySpace>
668
669
670 template <typename MemorySpace>
671 void
674 {
676 column_space_map =
677 Utilities::Trilinos::internal::make_rcp<MapType>(*sp.column_space_map);
678 graph = Utilities::Trilinos::internal::make_rcp<GraphType>(*sp.graph);
679
680 if (sp.nonlocal_graph.get() != nullptr)
681 nonlocal_graph = Utilities::Trilinos::internal::make_rcp<GraphType>(
682 *sp.nonlocal_graph);
683 else
684 nonlocal_graph.reset();
685 }
686
687
688
689 template <typename MemorySpace>
690 template <typename SparsityPatternType>
691 void
692 SparsityPattern<MemorySpace>::copy_from(const SparsityPatternType &sp)
693 {
694 SparsityPatternBase::resize(sp.n_rows(), sp.n_cols());
695 Teuchos::RCP<MapType> rows =
696 Utilities::Trilinos::internal::make_rcp<MapType>(
698 0,
700 Teuchos::RCP<MapType> columns =
701 Utilities::Trilinos::internal::make_rcp<MapType>(
703 0,
705
706 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
707 rows, columns, sp, false, column_space_map, graph, nonlocal_graph);
708 }
709
711
712 template <typename MemorySpace>
713 void
715 {
717 // When we clear the matrix, reset
718 // the pointer and generate an
719 // empty sparsity pattern.
720 column_space_map = Utilities::Trilinos::internal::make_rcp<MapType>(
724 graph =
725 Utilities::Trilinos::internal::make_rcp<GraphType>(column_space_map,
726 column_space_map,
727 0);
728 graph->fillComplete();
729
730 nonlocal_graph.reset();
732
733
734
735 template <typename MemorySpace>
736 void
738 {
739 Assert(column_space_map.get(), ExcInternalError());
740 if (nonlocal_graph.get() != nullptr)
741 {
742# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
743 if (nonlocal_graph->getRowMap()->getLocalNumElements() > 0 &&
744 column_space_map->getGlobalNumElements() > 0)
745# else
746 if (nonlocal_graph->getRowMap()->getNodeNumElements() > 0 &&
747 column_space_map->getGlobalNumElements() > 0)
748# endif
749 {
750 // Insert dummy element at (row, column) that corresponds to row 0
751 // in local index counting.
753 nonlocal_graph->getRowMap()->getGlobalElement(0);
755
756 // in case we have a square sparsity pattern, add the entry on the
757 // diagonal
758 if (column_space_map->getGlobalNumElements() ==
759 graph->getRangeMap()->getGlobalNumElements())
760 column = row;
761 // if not, take a column index that we have ourselves since we
762 // know for sure it is there (and it will not create spurious
763 // messages to many ranks like putting index 0 on many
764 // processors)
765# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
766 else if (column_space_map->getLocalNumElements() > 0)
767# else
768 else if (column_space_map->getNodeNumElements() > 0)
769# endif
770 column = column_space_map->getGlobalElement(0);
771 nonlocal_graph->insertGlobalIndices(row, 1, &column);
772 }
773# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
774 Assert(nonlocal_graph->getRowMap()->getLocalNumElements() == 0 ||
775 column_space_map->getGlobalNumElements() == 0,
777# else
778 Assert(nonlocal_graph->getRowMap()->getNodeNumElements() == 0 ||
779 column_space_map->getGlobalNumElements() == 0,
781# endif
782
783 nonlocal_graph->fillComplete(column_space_map, graph->getRowMap());
784 }
785 graph->fillComplete(column_space_map, graph->getRowMap());
786
787 // Check consistency between the sizes set at the beginning and what
788 // Trilinos stores:
789 using namespace deal_II_exceptions::internals;
790 AssertDimension(n_rows(), graph->getGlobalNumRows());
791 AssertDimension(n_cols(), graph->getGlobalNumCols());
792 }
793
794
795
796 template <typename MemorySpace>
797 bool
799 {
800 return graph->getRowMap()->getLocalElement(i) !=
801 Teuchos::OrdinalTraits<int>::invalid();
802 }
803
804
805
806 template <typename MemorySpace>
807 bool
809 const size_type j) const
810 {
811# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
812 if (!row_is_stored_locally(i))
813 return false;
814
815 // Extract local indices in the matrix.
816 const auto trilinos_i = graph->getRowMap()->getLocalElement(i);
817 const auto trilinos_j = graph->getColMap()->getLocalElement(j);
818
819 typename GraphType::local_inds_host_view_type col_indices;
820
821 // Generate the view.
822 graph->getLocalRowView(trilinos_i, col_indices);
823
824 // Search the index
825 const size_type local_col_index =
826 std::find(col_indices.data(),
827 col_indices.data() + col_indices.size(),
828 trilinos_j) -
829 col_indices.data();
830
831 return static_cast<size_t>(local_col_index) != col_indices.size();
832# else
833 (void)i;
834 (void)j;
836 return false;
837# endif
838 }
839
841
842 template <typename MemorySpace>
845 {
846 size_type local_b = 0;
847 for (int i = 0; i < static_cast<int>(local_size()); ++i)
848 {
849# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
850 typename GraphType::local_inds_host_view_type indices;
851# else
852 Teuchos::ArrayView<const int> indices;
853# endif
854
855 graph->getLocalRowView(i, indices);
856 const auto num_entries = indices.size();
857 for (unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
858 ++j)
859 {
860 if (static_cast<size_type>(std::abs(i - indices[j])) > local_b)
861 local_b = std::abs(i - indices[j]);
862 }
863 }
864
866 Utilities::MPI::max(local_b,
868 graph->getComm()));
869 return static_cast<size_type>(global_b);
870 }
871
872
873
874 template <typename MemorySpace>
875 unsigned int
877 {
878# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
879 return graph->getLocalNumRows();
880# else
881 return graph->getNodeNumRows();
882# endif
883 }
884
885
886
887 template <typename MemorySpace>
888 std::pair<typename SparsityPattern<MemorySpace>::size_type,
891 {
892 const size_type begin = graph->getRowMap()->getMinGlobalIndex();
893 const size_type end = graph->getRowMap()->getMaxGlobalIndex() + 1;
894
895 return {begin, end};
896 }
897
898
899
900 template <typename MemorySpace>
901 std::uint64_t
903 {
904 return graph->getGlobalNumEntries();
905 }
906
907
908
909 template <typename MemorySpace>
910 unsigned int
912 {
913# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
914 return graph->getLocalMaxNumRowEntries();
915# else
916 return graph->getNodeMaxNumRowEntries();
917# endif
918 }
919
920
921
922 template <typename MemorySpace>
925 {
926 Assert(row < (size_type)n_rows(), ExcInternalError());
927
928 // Get a representation of the where the present row is located on
929 // the current processor
931 graph->getRowMap()->getLocalElement(row);
932
933 // On the processor who owns this row, we'll have a non-negative
934 // value for `local_row` and can ask for the length of the row.
935 if (local_row >= 0)
936 return graph->getNumEntriesInLocalRow(local_row);
937 else
938 return static_cast<size_type>(-1);
939 }
940
941
942
943 template <typename MemorySpace>
944 void
946 const ::types::global_dof_index &row,
948 const bool indices_are_sorted)
949 {
950 add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
951 }
952
953
954
955 template <typename MemorySpace>
956 Teuchos::RCP<const typename SparsityPattern<MemorySpace>::MapType>
958 {
959 return graph->getDomainMap();
960 }
961
962
963
964 template <typename MemorySpace>
965 Teuchos::RCP<const typename SparsityPattern<MemorySpace>::MapType>
967 {
968 return graph->getRangeMap();
969 }
970
971
972
973 template <typename MemorySpace>
976 {
978 graph->getRangeMap()->getComm());
979 }
980
981
982
983 template <typename MemorySpace>
984 Teuchos::RCP<const Teuchos::Comm<int>>
987 return graph->getRangeMap()->getComm();
988 }
989
990
991
992 // As of now, no particularly neat
993 // output is generated in case of
994 // multiple processors.
995 template <typename MemorySpace>
996 void
998 std::ostream &out,
999 const bool write_extended_trilinos_info) const
1000 {
1001 if (write_extended_trilinos_info)
1002 out << *graph;
1003 else
1004 {
1005# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
1006 for (unsigned int i = 0; i < graph->getLocalNumRows(); ++i)
1007# else
1008 for (unsigned int i = 0; i < graph->getNodeNumRows(); ++i)
1009# endif
1010 {
1011# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1012 typename GraphType::local_inds_host_view_type indices;
1013# else
1014 Teuchos::ArrayView<const int> indices;
1015# endif
1016 graph->getLocalRowView(i, indices);
1017 int num_entries = indices.size();
1018 for (int j = 0; j < num_entries; ++j)
1019 out << "(" << graph->getRowMap()->getGlobalElement(i) << ","
1020 << graph->getColMap()->getGlobalElement(indices[j]) << ") "
1021 << std::endl;
1022 }
1023 }
1024
1025 AssertThrow(out.fail() == false, ExcIO());
1026 }
1027
1028
1029
1030 template <typename MemorySpace>
1031 void
1033 {
1034 Assert(graph->isFillComplete() == true, ExcInternalError());
1035
1036 for (unsigned int row = 0; row < local_size(); ++row)
1037 {
1038# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1039 typename GraphType::local_inds_host_view_type indices;
1040# else
1041 Teuchos::ArrayView<const int> indices;
1042# endif
1043 graph->getLocalRowView(row, indices);
1044 int num_entries = indices.size();
1045
1046 Assert(num_entries >= 0, ExcInternalError());
1047 // avoid sign comparison warning
1048 const ::types::signed_global_dof_index num_entries_ =
1049 num_entries;
1050 for (::types::signed_global_dof_index j = 0; j < num_entries_;
1051 ++j)
1052 // while matrix entries are usually
1053 // written (i,j), with i vertical and
1054 // j horizontal, gnuplot output is
1055 // x-y, that is we have to exchange
1056 // the order of output
1057 out << static_cast<int>(
1058 graph->getColMap()->getGlobalElement(indices[j]))
1059 << " "
1060 << -static_cast<int>(graph->getRowMap()->getGlobalElement(row))
1061 << std::endl;
1062 }
1063
1064 AssertThrow(out.fail() == false, ExcIO());
1065 }
1066
1067 // TODO: Implement!
1068 template <typename MemorySpace>
1069 std::size_t
1075
1076
1077# ifndef DOXYGEN
1078 // explicit instantiations
1080
1081 template void
1083 const ::SparsityPattern &);
1084 template void
1086 const ::DynamicSparsityPattern &);
1087
1088 template void
1090 const IndexSet &,
1091 const ::SparsityPattern &,
1092 const MPI_Comm,
1093 bool);
1094 template void
1096 const IndexSet &,
1097 const ::DynamicSparsityPattern &,
1098 const MPI_Comm,
1099 bool);
1100
1101
1102 template void
1104 const IndexSet &,
1105 const IndexSet &,
1106 const ::SparsityPattern &,
1107 const MPI_Comm,
1108 bool);
1109 template void
1111 const IndexSet &,
1112 const IndexSet &,
1113 const ::DynamicSparsityPattern &,
1114 const MPI_Comm,
1115 bool);
1116# endif
1117
1118 } // namespace TpetraWrappers
1119
1120} // namespace LinearAlgebra
1121
1123
1124#endif // DEAL_II_TRILINOS_WITH_TPETRA
iterator begin() const
Definition array_view.h:702
iterator end() const
Definition array_view.h:711
size_type size() const
Definition index_set.h:1765
size_type n_elements() const
Definition index_set.h:1923
void subtract_set(const IndexSet &other)
Definition index_set.cc:470
Teuchos::RCP< Tpetra::Map< int, types::signed_global_dof_index > > make_tpetra_map_rcp(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:962
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
Teuchos::RCP< const Teuchos::Comm< int > > get_teuchos_mpi_communicator() 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
Tpetra::CrsGraph< int, ::types::signed_global_dof_index, NodeType > GraphType
bool exists(const size_type i, const size_type j) const
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
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
Tpetra::Map< int, ::types::signed_global_dof_index, NodeType > MapType
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
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:1193
Tpetra::CrsGraph< int, ::types::signed_global_dof_index, typename SparsityPattern< MemorySpace >::NodeType > GraphType
Tpetra::Map< int, ::types::signed_global_dof_index, typename SparsityPattern< MemorySpace >::NodeType > MapType
void reinit_sp(const Teuchos::RCP< MapType< MemorySpace > > &row_map, const Teuchos::RCP< MapType< MemorySpace > > &col_map, const size_type< MemorySpace > n_entries_per_row, Teuchos::RCP< MapType< MemorySpace > > &column_space_map, Teuchos::RCP< GraphType< MemorySpace > > &graph, Teuchos::RCP< GraphType< MemorySpace > > &nonlocal_graph)
typename SparsityPattern< MemorySpace >::size_type size_type
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:92
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 > &)
int signed_global_dof_index
Definition types.h:92