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