Reference documentation for deal.II version GIT relicensing-677-g613a5db665 2024-05-19 13:30: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_sparse_matrix.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17
18#ifdef DEAL_II_WITH_TRILINOS
19
22
29
30# include <boost/container/small_vector.hpp>
31
32# ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT
33# include <EpetraExt_MatrixMatrix.h>
34# endif
35# include <Epetra_Export.h>
36# include <Teuchos_RCP.hpp>
37# include <ml_epetra_utils.h>
38# include <ml_struct.h>
39
40# include <memory>
41
43
44namespace TrilinosWrappers
45{
46 namespace internal
47 {
48 template <typename VectorType>
49 typename VectorType::value_type *
50 begin(VectorType &V)
51 {
52 return V.begin();
53 }
54
55 template <typename VectorType>
56 const typename VectorType::value_type *
57 begin(const VectorType &V)
58 {
59 return V.begin();
60 }
61
62 template <typename VectorType>
63 typename VectorType::value_type *
64 end(VectorType &V)
65 {
66 return V.end();
67 }
68
69 template <typename VectorType>
70 const typename VectorType::value_type *
71 end(const VectorType &V)
72 {
73 return V.end();
74 }
75
76 template <>
77 double *
82
83 template <>
84 const double *
86 {
87 return V.trilinos_vector()[0];
88 }
89
90 template <>
91 double *
93 {
94 return V.trilinos_vector()[0] + V.trilinos_vector().MyLength();
95 }
96
97 template <>
98 const double *
100 {
101 return V.trilinos_vector()[0] + V.trilinos_vector().MyLength();
102 }
103
104# ifdef DEAL_II_TRILINOS_WITH_TPETRA
105 template <typename Number>
106 Number *
108 {
109 return V.trilinos_vector().getDataNonConst().get();
110 }
111
112 template <typename Number>
113 const Number *
115 {
116 return V.trilinos_vector().getData().get();
117 }
118
119 template <typename Number>
120 Number *
122 {
123 return V.trilinos_vector().getDataNonConst().get() +
124 V.trilinos_vector().getLocalLength();
125 }
126
127 template <typename Number>
128 const Number *
130 {
131 return V.trilinos_vector().getData().get() +
132 V.trilinos_vector().getLocalLength();
133 }
134# endif
135 } // namespace internal
136
137
138 namespace SparseMatrixIterators
139 {
140 void
142 {
143 // if we are asked to visit the past-the-end line, then simply
144 // release all our caches and go on with life.
145 //
146 // do the same if the row we're supposed to visit is not locally
147 // owned. this is simply going to make non-locally owned rows
148 // look like they're empty
149 if ((this->a_row == matrix->m()) ||
150 (matrix->in_local_range(this->a_row) == false))
151 {
152 colnum_cache.reset();
153 value_cache.reset();
154
155 return;
156 }
157
158 // get a representation of the present row
159 int ncols;
161 matrix->row_length(this->a_row);
162 if (value_cache.get() == nullptr)
163 {
164 value_cache = std::make_shared<std::vector<TrilinosScalar>>(colnums);
165 colnum_cache = std::make_shared<std::vector<size_type>>(colnums);
166 }
167 else
168 {
169 value_cache->resize(colnums);
170 colnum_cache->resize(colnums);
171 }
172
173 const int ierr = matrix->trilinos_matrix().ExtractGlobalRowCopy(
174 this->a_row,
175 colnums,
176 ncols,
177 value_cache->data(),
178 reinterpret_cast<TrilinosWrappers::types::int_type *>(
179 colnum_cache->data()));
180 AssertDimension(ncols, colnums);
181 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
182
183 // copy it into our caches if the
184 // line isn't empty. if it is, then
185 // we've done something wrong, since
186 // we shouldn't have initialized an
187 // iterator for an empty line (what
188 // would it point to?)
189 }
190 } // namespace SparseMatrixIterators
191
192
193 // The constructor is actually the
194 // only point where we have to check
195 // whether we build a serial or a
196 // parallel Trilinos matrix.
197 // Actually, it does not even matter
198 // how many threads there are, but
199 // only if we use an MPI compiler or
200 // a standard compiler. So, even one
201 // thread on a configuration with
202 // MPI will still get a parallel
203 // interface.
205 : column_space_map(new Epetra_Map(0, 0, Utilities::Trilinos::comm_self()))
206 , matrix(
207 new Epetra_FECrsMatrix(View, *column_space_map, *column_space_map, 0))
208 , last_action(Zero)
209 , compressed(true)
210 {
211 matrix->FillComplete();
212 }
213
214
215
217 const size_type n,
218 const unsigned int n_max_entries_per_row)
219 : column_space_map(
220 new Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(n),
221 0,
222 Utilities::Trilinos::comm_self()))
223 ,
224
225 // on one processor only, we know how the
226 // columns of the matrix will be
227 // distributed (everything on one
228 // processor), so we can hand in this
229 // information to the constructor. we
230 // can't do so in parallel, where the
231 // information from columns is only
232 // available when entries have been added
233 matrix(new Epetra_FECrsMatrix(
234 Copy,
235 Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(m),
236 0,
237 Utilities::Trilinos::comm_self()),
238 *column_space_map,
239 n_max_entries_per_row,
240 false))
241 , last_action(Zero)
242 , compressed(false)
243 {}
244
245
246
248 const size_type n,
249 const std::vector<unsigned int> &n_entries_per_row)
250 : column_space_map(
251 new Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(n),
252 0,
253 Utilities::Trilinos::comm_self()))
254 , matrix(new Epetra_FECrsMatrix(
255 Copy,
256 Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(m),
257 0,
258 Utilities::Trilinos::comm_self()),
259 *column_space_map,
260 reinterpret_cast<int *>(
261 const_cast<unsigned int *>(n_entries_per_row.data())),
262 false))
263 , last_action(Zero)
264 , compressed(false)
265 {}
266
267
268
269 SparseMatrix::SparseMatrix(const IndexSet &parallel_partitioning,
270 const MPI_Comm communicator,
271 const unsigned int n_max_entries_per_row)
272 : column_space_map(new Epetra_Map(
273 parallel_partitioning.make_trilinos_map(communicator, false)))
274 , matrix(new Epetra_FECrsMatrix(Copy,
275 *column_space_map,
276 n_max_entries_per_row,
277 false))
278 , last_action(Zero)
279 , compressed(false)
280 {}
281
282
283
284 SparseMatrix::SparseMatrix(const IndexSet &parallel_partitioning,
285 const MPI_Comm communicator,
286 const std::vector<unsigned int> &n_entries_per_row)
287 : column_space_map(new Epetra_Map(
288 parallel_partitioning.make_trilinos_map(communicator, false)))
289 , matrix(new Epetra_FECrsMatrix(Copy,
290 *column_space_map,
291 reinterpret_cast<int *>(
292 const_cast<unsigned int *>(
293 n_entries_per_row.data())),
294 false))
295 , last_action(Zero)
296 , compressed(false)
297 {}
298
299
300
301 SparseMatrix::SparseMatrix(const IndexSet &row_parallel_partitioning,
302 const IndexSet &col_parallel_partitioning,
303 const MPI_Comm communicator,
304 const size_type n_max_entries_per_row)
305 : column_space_map(new Epetra_Map(
306 col_parallel_partitioning.make_trilinos_map(communicator, false)))
307 , matrix(new Epetra_FECrsMatrix(
308 Copy,
309 row_parallel_partitioning.make_trilinos_map(communicator, false),
310 n_max_entries_per_row,
311 false))
312 , last_action(Zero)
313 , compressed(false)
314 {}
315
316
317
318 SparseMatrix::SparseMatrix(const IndexSet &row_parallel_partitioning,
319 const IndexSet &col_parallel_partitioning,
320 const MPI_Comm communicator,
321 const std::vector<unsigned int> &n_entries_per_row)
322 : column_space_map(new Epetra_Map(
323 col_parallel_partitioning.make_trilinos_map(communicator, false)))
324 , matrix(new Epetra_FECrsMatrix(
325 Copy,
326 row_parallel_partitioning.make_trilinos_map(communicator, false),
327 reinterpret_cast<int *>(
328 const_cast<unsigned int *>(n_entries_per_row.data())),
329 false))
330 , last_action(Zero)
331 , compressed(false)
332 {}
333
334
335
337 : column_space_map(new Epetra_Map(sparsity_pattern.domain_partitioner()))
338 , matrix(
339 new Epetra_FECrsMatrix(Copy,
340 sparsity_pattern.trilinos_sparsity_pattern(),
341 false))
342 , last_action(Zero)
343 , compressed(true)
344 {
345 Assert(sparsity_pattern.trilinos_sparsity_pattern().Filled() == true,
347 "The Trilinos sparsity pattern has not been compressed."));
349 }
350
351
352
354 : column_space_map(std::move(other.column_space_map))
355 , matrix(std::move(other.matrix))
356 , nonlocal_matrix(std::move(other.nonlocal_matrix))
357 , nonlocal_matrix_exporter(std::move(other.nonlocal_matrix_exporter))
358 , last_action(other.last_action)
359 , compressed(other.compressed)
360 {
361 other.last_action = Zero;
362 other.compressed = false;
363 }
364
365
366
367 void
369 {
370 if (this == &rhs)
371 return;
372
373 nonlocal_matrix.reset();
375
376 // check whether we need to update the whole matrix layout (we have
377 // different maps or if we detect a row where the columns of the two
378 // matrices do not match)
379 bool needs_deep_copy =
380 !matrix->RowMap().SameAs(rhs.matrix->RowMap()) ||
381 !matrix->ColMap().SameAs(rhs.matrix->ColMap()) ||
382 !matrix->DomainMap().SameAs(rhs.matrix->DomainMap()) ||
384 if (!needs_deep_copy)
385 {
386 // Try to copy all the rows of the matrix one by one. In case of error
387 // (i.e., the column indices are different), we need to abort and blow
388 // away the matrix.
389 for (const auto row : locally_owned_range_indices())
390 {
391 const int row_local = matrix->RowMap().LID(
392 static_cast<TrilinosWrappers::types::int_type>(row));
393 Assert((row_local >= 0), ExcAccessToNonlocalRow(row));
394
395 int n_entries, rhs_n_entries;
396 TrilinosScalar *value_ptr, *rhs_value_ptr;
397 int *index_ptr, *rhs_index_ptr;
398 int ierr = rhs.matrix->ExtractMyRowView(row_local,
399 rhs_n_entries,
400 rhs_value_ptr,
401 rhs_index_ptr);
402 (void)ierr;
403 Assert(ierr == 0, ExcTrilinosError(ierr));
404
405 ierr = matrix->ExtractMyRowView(row_local,
406 n_entries,
407 value_ptr,
408 index_ptr);
409 Assert(ierr == 0, ExcTrilinosError(ierr));
410
411 if (n_entries != rhs_n_entries ||
412 std::memcmp(static_cast<void *>(index_ptr),
413 static_cast<void *>(rhs_index_ptr),
414 sizeof(int) * n_entries) != 0)
415 {
416 needs_deep_copy = true;
417 break;
418 }
419
420 for (int i = 0; i < n_entries; ++i)
421 value_ptr[i] = rhs_value_ptr[i];
422 }
423 }
424
425 if (needs_deep_copy)
426 {
428 std::make_unique<Epetra_Map>(rhs.trilinos_matrix().DomainMap());
429
430 // release memory before reallocation
431 matrix = std::make_unique<Epetra_FECrsMatrix>(*rhs.matrix);
432
433 matrix->FillComplete(*column_space_map, matrix->RowMap());
434 }
435
436 if (rhs.nonlocal_matrix.get() != nullptr)
438 std::make_unique<Epetra_CrsMatrix>(Copy, rhs.nonlocal_matrix->Graph());
439 }
440
441
442
443 namespace
444 {
445 template <typename SparsityPatternType>
446 void
447 reinit_matrix(const IndexSet &row_parallel_partitioning,
448 const IndexSet &column_parallel_partitioning,
449 const SparsityPatternType &sparsity_pattern,
450 const bool exchange_data,
451 const MPI_Comm communicator,
452 std::unique_ptr<Epetra_Map> &column_space_map,
453 std::unique_ptr<Epetra_FECrsMatrix> &matrix,
454 std::unique_ptr<Epetra_CrsMatrix> &nonlocal_matrix,
455 std::unique_ptr<Epetra_Export> &nonlocal_matrix_exporter)
456 {
457 // release memory before reallocation
458 matrix.reset();
459 nonlocal_matrix.reset();
460 nonlocal_matrix_exporter.reset();
461
462 column_space_map = std::make_unique<Epetra_Map>(
463 column_parallel_partitioning.make_trilinos_map(communicator, false));
464
465 if (column_space_map->Comm().MyPID() == 0)
466 {
467 AssertDimension(sparsity_pattern.n_rows(),
468 row_parallel_partitioning.size());
469 AssertDimension(sparsity_pattern.n_cols(),
470 column_parallel_partitioning.size());
471 }
472
473 Epetra_Map row_space_map =
474 row_parallel_partitioning.make_trilinos_map(communicator, false);
475
476 // if we want to exchange data, build a usual Trilinos sparsity pattern
477 // and let that handle the exchange. otherwise, manually create a
478 // CrsGraph, which consumes considerably less memory because it can set
479 // correct number of indices right from the start
480 if (exchange_data)
481 {
482 SparsityPattern trilinos_sparsity;
483 trilinos_sparsity.reinit(row_parallel_partitioning,
484 column_parallel_partitioning,
485 sparsity_pattern,
486 communicator,
487 exchange_data);
488 matrix = std::make_unique<Epetra_FECrsMatrix>(
489 Copy, trilinos_sparsity.trilinos_sparsity_pattern(), false);
490
491 return;
492 }
493
495 row_space_map),
497 row_space_map) +
498 1;
499 std::vector<int> n_entries_per_row(last_row - first_row);
500
501 for (SparseMatrix::size_type row = first_row; row < last_row; ++row)
502 n_entries_per_row[row - first_row] = sparsity_pattern.row_length(row);
503
504 // The deal.II notation of a Sparsity pattern corresponds to the Epetra
505 // concept of a Graph. Hence, we generate a graph by copying the
506 // sparsity pattern into it, and then build up the matrix from the
507 // graph. This is considerable faster than directly filling elements
508 // into the matrix. Moreover, it consumes less memory, since the
509 // internal reordering is done on ints only, and we can leave the
510 // doubles aside.
511
512 // for more than one processor, need to specify only row map first and
513 // let the matrix entries decide about the column map (which says which
514 // columns are present in the matrix, not to be confused with the
515 // col_map that tells how the domain dofs of the matrix will be
516 // distributed). for only one processor, we can directly assign the
517 // columns as well. Compare this with bug # 4123 in the Sandia Bugzilla.
518 std::unique_ptr<Epetra_CrsGraph> graph;
519 if (row_space_map.Comm().NumProc() > 1)
520 graph = std::make_unique<Epetra_CrsGraph>(Copy,
521 row_space_map,
522 n_entries_per_row.data(),
523 true);
524 else
525 graph = std::make_unique<Epetra_CrsGraph>(Copy,
526 row_space_map,
527 *column_space_map,
528 n_entries_per_row.data(),
529 true);
530
531 // This functions assumes that the sparsity pattern sits on all
532 // processors (completely). The parallel version uses an Epetra graph
533 // that is already distributed.
534
535 // now insert the indices
536 std::vector<TrilinosWrappers::types::int_type> row_indices;
537
538 for (SparseMatrix::size_type row = first_row; row < last_row; ++row)
539 {
540 const int row_length = sparsity_pattern.row_length(row);
541 if (row_length == 0)
542 continue;
543
544 row_indices.resize(row_length, -1);
545 {
546 typename SparsityPatternType::iterator p =
547 sparsity_pattern.begin(row);
548 for (SparseMatrix::size_type col = 0;
549 p != sparsity_pattern.end(row);
550 ++p, ++col)
551 row_indices[col] = p->column();
552 }
553 graph->Epetra_CrsGraph::InsertGlobalIndices(row,
554 row_length,
555 row_indices.data());
556 }
557
558 // Eventually, optimize the graph structure (sort indices, make memory
559 // contiguous, etc). note that the documentation of the function indeed
560 // states that we first need to provide the column (domain) map and then
561 // the row (range) map
562 graph->FillComplete(*column_space_map, row_space_map);
563 graph->OptimizeStorage();
564
565 // check whether we got the number of columns right.
566 AssertDimension(sparsity_pattern.n_cols(),
568 (void)n_global_cols;
569
570 // And now finally generate the matrix.
571 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph, false);
572 }
573
574
575
576 // for the non-local graph, we need to circumvent the problem that some
577 // processors will not add into the non-local graph at all: We do not want
578 // to insert dummy elements on >5000 processors because that gets very
579 // slow. Thus, we set a flag in Epetra_CrsGraph that sets the correct
580 // flag. Since it is protected, we need to expose this information by
581 // deriving a class from Epetra_CrsGraph for the purpose of creating the
582 // data structure
583 class Epetra_CrsGraphMod : public Epetra_CrsGraph
584 {
585 public:
586 Epetra_CrsGraphMod(const Epetra_Map &row_map,
587 const int *n_entries_per_row)
588 : Epetra_CrsGraph(Copy, row_map, n_entries_per_row, true)
589 {}
590
591 void
592 SetIndicesAreGlobal()
593 {
594 this->Epetra_CrsGraph::SetIndicesAreGlobal(true);
595 }
596 };
597
598
599
600 // specialization for DynamicSparsityPattern which can provide us with
601 // more information about the non-locally owned rows
602 template <>
603 void
604 reinit_matrix(const IndexSet &row_parallel_partitioning,
605 const IndexSet &column_parallel_partitioning,
606 const DynamicSparsityPattern &sparsity_pattern,
607 const bool exchange_data,
608 const MPI_Comm communicator,
609 std::unique_ptr<Epetra_Map> &column_space_map,
610 std::unique_ptr<Epetra_FECrsMatrix> &matrix,
611 std::unique_ptr<Epetra_CrsMatrix> &nonlocal_matrix,
612 std::unique_ptr<Epetra_Export> &nonlocal_matrix_exporter)
613 {
614 matrix.reset();
615 nonlocal_matrix.reset();
616 nonlocal_matrix_exporter.reset();
617
618 column_space_map = std::make_unique<Epetra_Map>(
619 column_parallel_partitioning.make_trilinos_map(communicator, false));
620
621 AssertDimension(sparsity_pattern.n_rows(),
622 row_parallel_partitioning.size());
623 AssertDimension(sparsity_pattern.n_cols(),
624 column_parallel_partitioning.size());
625
626 Epetra_Map row_space_map =
627 row_parallel_partitioning.make_trilinos_map(communicator, false);
628
629 IndexSet relevant_rows(sparsity_pattern.row_index_set());
630 // serial case
631 if (relevant_rows.size() == 0)
632 {
633 relevant_rows.set_size(
635 relevant_rows.add_range(
636 0, TrilinosWrappers::n_global_elements(row_space_map));
637 }
638 relevant_rows.compress();
639 Assert(relevant_rows.n_elements() >=
640 static_cast<unsigned int>(row_space_map.NumMyElements()),
642 "Locally relevant rows of sparsity pattern must contain "
643 "all locally owned rows"));
644
645 // check whether the relevant rows correspond to exactly the same map as
646 // the owned rows. In that case, do not create the nonlocal graph and
647 // fill the columns by demand
648 const bool have_ghost_rows = [&]() {
649 const std::vector<::types::global_dof_index> indices =
650 relevant_rows.get_index_vector();
651 Epetra_Map relevant_map(
653 TrilinosWrappers::types::int_type(relevant_rows.n_elements()),
654 (indices.empty() ?
655 nullptr :
656 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
657 indices.data())),
658 0,
659 row_space_map.Comm());
660 return !relevant_map.SameAs(row_space_map);
661 }();
662
663 std::vector<TrilinosWrappers::types::int_type> ghost_rows;
664 std::vector<int> n_entries_per_row(row_space_map.NumMyElements());
665 std::vector<int> n_entries_per_ghost_row;
666 {
668 for (const auto global_row : relevant_rows)
669 {
670 if (row_space_map.MyGID(
671 static_cast<TrilinosWrappers::types::int_type>(global_row)))
672 n_entries_per_row[own++] =
673 sparsity_pattern.row_length(global_row);
674 else if (sparsity_pattern.row_length(global_row) > 0)
675 {
676 ghost_rows.push_back(global_row);
677 n_entries_per_ghost_row.push_back(
678 sparsity_pattern.row_length(global_row));
679 }
680 }
681 }
682
683 Epetra_Map off_processor_map(-1,
684 ghost_rows.size(),
685 (ghost_rows.size() > 0) ?
686 (ghost_rows.data()) :
687 nullptr,
688 0,
689 row_space_map.Comm());
690
691 std::unique_ptr<Epetra_CrsGraph> graph;
692 std::unique_ptr<Epetra_CrsGraphMod> nonlocal_graph;
693 if (row_space_map.Comm().NumProc() > 1)
694 {
695 graph =
696 std::make_unique<Epetra_CrsGraph>(Copy,
697 row_space_map,
698 (n_entries_per_row.size() > 0) ?
699 (n_entries_per_row.data()) :
700 nullptr,
701 exchange_data ? false : true);
702 if (have_ghost_rows == true)
703 nonlocal_graph = std::make_unique<Epetra_CrsGraphMod>(
704 off_processor_map, n_entries_per_ghost_row.data());
705 }
706 else
707 graph =
708 std::make_unique<Epetra_CrsGraph>(Copy,
709 row_space_map,
710 *column_space_map,
711 (n_entries_per_row.size() > 0) ?
712 (n_entries_per_row.data()) :
713 nullptr,
714 true);
715
716 // now insert the indices, select between the right matrix
717 std::vector<TrilinosWrappers::types::int_type> row_indices;
718
719 for (const auto global_row : relevant_rows)
720 {
721 const int row_length = sparsity_pattern.row_length(global_row);
722 if (row_length == 0)
723 continue;
724
725 row_indices.resize(row_length, -1);
726 for (int col = 0; col < row_length; ++col)
727 row_indices[col] = sparsity_pattern.column_number(global_row, col);
728
729 if (row_space_map.MyGID(
730 static_cast<TrilinosWrappers::types::int_type>(global_row)))
731 graph->InsertGlobalIndices(global_row,
732 row_length,
733 row_indices.data());
734 else
735 {
736 Assert(nonlocal_graph.get() != nullptr, ExcInternalError());
737 nonlocal_graph->InsertGlobalIndices(global_row,
738 row_length,
739 row_indices.data());
740 }
741 }
742
743 // finalize nonlocal graph and create nonlocal matrix
744 if (nonlocal_graph.get() != nullptr)
745 {
746 // must make sure the IndicesAreGlobal flag is set on all processors
747 // because some processors might not call InsertGlobalIndices (and
748 // we do not want to insert dummy indices on all processors for
749 // large-scale simulations due to the bad impact on performance)
750 nonlocal_graph->SetIndicesAreGlobal();
751 Assert(nonlocal_graph->IndicesAreGlobal() == true,
753 nonlocal_graph->FillComplete(*column_space_map, row_space_map);
754 nonlocal_graph->OptimizeStorage();
755
756 // insert data from nonlocal graph into the final sparsity pattern
757 if (exchange_data)
758 {
759 Epetra_Export exporter(nonlocal_graph->RowMap(), row_space_map);
760 int ierr = graph->Export(*nonlocal_graph, exporter, Add);
761 (void)ierr;
762 Assert(ierr == 0, ExcTrilinosError(ierr));
763 }
764
765 nonlocal_matrix =
766 std::make_unique<Epetra_CrsMatrix>(Copy, *nonlocal_graph);
767 }
768
769 graph->FillComplete(*column_space_map, row_space_map);
770 graph->OptimizeStorage();
771
772 AssertDimension(sparsity_pattern.n_cols(),
774
775 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph, false);
776 }
777 } // namespace
778
779
780
781 template <typename SparsityPatternType>
782 void
783 SparseMatrix::reinit(const SparsityPatternType &sparsity_pattern)
784 {
785 reinit_matrix(complete_index_set(sparsity_pattern.n_rows()),
786 complete_index_set(sparsity_pattern.n_cols()),
787 sparsity_pattern,
788 false,
789 MPI_COMM_SELF,
791 matrix,
794 }
795
796
797
798 template <typename SparsityPatternType>
799 inline std::enable_if_t<
800 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
801 SparseMatrix::reinit(const IndexSet &row_parallel_partitioning,
802 const IndexSet &col_parallel_partitioning,
803 const SparsityPatternType &sparsity_pattern,
804 const MPI_Comm communicator,
805 const bool exchange_data)
806 {
807 reinit_matrix(row_parallel_partitioning,
808 col_parallel_partitioning,
809 sparsity_pattern,
810 exchange_data,
811 communicator,
813 matrix,
816
817 // In the end, the matrix needs to be compressed in order to be really
818 // ready.
819 last_action = Zero;
821 }
822
823
824
825 void
826 SparseMatrix::reinit(const SparsityPattern &sparsity_pattern)
827 {
828 matrix.reset();
830
831 // reinit with a (parallel) Trilinos sparsity pattern.
833 std::make_unique<Epetra_Map>(sparsity_pattern.domain_partitioner());
834 matrix = std::make_unique<Epetra_FECrsMatrix>(
835 Copy, sparsity_pattern.trilinos_sparsity_pattern(), false);
836
837 if (sparsity_pattern.nonlocal_graph.get() != nullptr)
839 std::make_unique<Epetra_CrsMatrix>(Copy,
840 *sparsity_pattern.nonlocal_graph);
841 else
842 nonlocal_matrix.reset();
843
844 last_action = Zero;
846 }
847
848
849
850 void
851 SparseMatrix::reinit(const SparseMatrix &sparse_matrix)
852 {
853 if (this == &sparse_matrix)
854 return;
855
857 std::make_unique<Epetra_Map>(sparse_matrix.trilinos_matrix().DomainMap());
858 matrix.reset();
860 matrix = std::make_unique<Epetra_FECrsMatrix>(
861 Copy, sparse_matrix.trilinos_sparsity_pattern(), false);
862
863 if (sparse_matrix.nonlocal_matrix != nullptr)
864 nonlocal_matrix = std::make_unique<Epetra_CrsMatrix>(
865 Copy, sparse_matrix.nonlocal_matrix->Graph());
866 else
867 nonlocal_matrix.reset();
868
869 last_action = Zero;
871 }
872
873
874
875 template <typename number>
876 inline void
878 const IndexSet &row_parallel_partitioning,
879 const IndexSet &col_parallel_partitioning,
880 const ::SparseMatrix<number> &dealii_sparse_matrix,
881 const MPI_Comm communicator,
882 const double drop_tolerance,
883 const bool copy_values,
884 const ::SparsityPattern *use_this_sparsity)
885 {
886 if (copy_values == false)
887 {
888 // in case we do not copy values, just
889 // call the other function.
890 if (use_this_sparsity == nullptr)
891 reinit(row_parallel_partitioning,
892 col_parallel_partitioning,
893 dealii_sparse_matrix.get_sparsity_pattern(),
894 communicator,
895 false);
896 else
897 reinit(row_parallel_partitioning,
898 col_parallel_partitioning,
899 *use_this_sparsity,
900 communicator,
901 false);
902 return;
903 }
904
905 const size_type n_rows = dealii_sparse_matrix.m();
906
907 AssertDimension(row_parallel_partitioning.size(), n_rows);
908 AssertDimension(col_parallel_partitioning.size(), dealii_sparse_matrix.n());
909
910 const ::SparsityPattern &sparsity_pattern =
911 (use_this_sparsity != nullptr) ?
912 *use_this_sparsity :
913 dealii_sparse_matrix.get_sparsity_pattern();
914
915 if (matrix.get() == nullptr || m() != n_rows ||
916 n_nonzero_elements() != sparsity_pattern.n_nonzero_elements())
917 {
918 reinit(row_parallel_partitioning,
919 col_parallel_partitioning,
920 sparsity_pattern,
921 communicator,
922 false);
923 }
924
925 // fill the values. the same as above: go through all rows of the
926 // matrix, and then all columns. since the sparsity patterns of the
927 // input matrix and the specified sparsity pattern might be different,
928 // need to go through the row for both these sparsity structures
929 // simultaneously in order to really set the correct values.
930 size_type maximum_row_length = matrix->MaxNumEntries();
931 std::vector<size_type> row_indices(maximum_row_length);
932 std::vector<TrilinosScalar> values(maximum_row_length);
933
934 for (size_type row = 0; row < n_rows; ++row)
935 // see if the row is locally stored on this processor
936 if (row_parallel_partitioning.is_element(row) == true)
937 {
938 ::SparsityPattern::iterator select_index =
939 sparsity_pattern.begin(row);
940 typename ::SparseMatrix<number>::const_iterator it =
941 dealii_sparse_matrix.begin(row);
942 size_type col = 0;
943 if (sparsity_pattern.n_rows() == sparsity_pattern.n_cols())
944 {
945 // optimized diagonal
946 AssertDimension(it->column(), row);
947 if (std::fabs(it->value()) > drop_tolerance)
948 {
949 values[col] = it->value();
950 row_indices[col++] = it->column();
951 }
952 ++select_index;
953 ++it;
954 }
955
956 while (it != dealii_sparse_matrix.end(row) &&
957 select_index != sparsity_pattern.end(row))
958 {
959 while (select_index->column() < it->column() &&
960 select_index != sparsity_pattern.end(row))
961 ++select_index;
962 while (it->column() < select_index->column() &&
963 it != dealii_sparse_matrix.end(row))
964 ++it;
965
966 if (it == dealii_sparse_matrix.end(row))
967 break;
968 if (std::fabs(it->value()) > drop_tolerance)
969 {
970 values[col] = it->value();
971 row_indices[col++] = it->column();
972 }
973 ++select_index;
974 ++it;
975 }
976 set(row,
977 col,
978 reinterpret_cast<size_type *>(row_indices.data()),
979 values.data(),
980 false);
981 }
983 }
984
985
986
987 template <typename number>
988 void
990 const ::SparseMatrix<number> &dealii_sparse_matrix,
991 const double drop_tolerance,
992 const bool copy_values,
993 const ::SparsityPattern *use_this_sparsity)
994 {
995 reinit(complete_index_set(dealii_sparse_matrix.m()),
996 complete_index_set(dealii_sparse_matrix.n()),
997 dealii_sparse_matrix,
998 MPI_COMM_SELF,
999 drop_tolerance,
1000 copy_values,
1001 use_this_sparsity);
1002 }
1003
1004
1005
1006 void
1007 SparseMatrix::reinit(const Epetra_CrsMatrix &input_matrix,
1008 const bool copy_values)
1009 {
1010 Assert(input_matrix.Filled() == true,
1011 ExcMessage("Input CrsMatrix has not called FillComplete()!"));
1012
1013 column_space_map = std::make_unique<Epetra_Map>(input_matrix.DomainMap());
1014
1015 const Epetra_CrsGraph *graph = &input_matrix.Graph();
1016
1017 nonlocal_matrix.reset();
1019 matrix.reset();
1020 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph, false);
1021
1022 matrix->FillComplete(*column_space_map, input_matrix.RangeMap(), true);
1023
1024 if (copy_values == true)
1025 {
1026 // point to the first data entry in the two
1027 // matrices and copy the content
1028 const TrilinosScalar *in_values = input_matrix[0];
1029 TrilinosScalar *values = (*matrix)[0];
1030 const size_type my_nonzeros = input_matrix.NumMyNonzeros();
1031 std::memcpy(values, in_values, my_nonzeros * sizeof(TrilinosScalar));
1032 }
1033
1034 last_action = Zero;
1036 }
1037
1038
1039
1040 void
1042 {
1043 Epetra_CombineMode mode = last_action;
1044 if (last_action == Zero)
1045 {
1046 if ((operation == VectorOperation::add) ||
1047 (operation == VectorOperation::unknown))
1048 mode = Add;
1049 else if (operation == VectorOperation::insert)
1050 mode = Insert;
1051 else
1052 Assert(
1053 false,
1054 ExcMessage(
1055 "compress() can only be called with VectorOperation add, insert, or unknown"));
1056 }
1057 else
1058 {
1059 Assert(
1060 ((last_action == Add) && (operation != VectorOperation::insert)) ||
1061 ((last_action == Insert) && (operation != VectorOperation::add)),
1062 ExcMessage("Operation and argument to compress() do not match"));
1063 }
1064
1065 // flush buffers
1066 int ierr;
1067 if (nonlocal_matrix.get() != nullptr && mode == Add)
1068 {
1069 // do only export in case of an add() operation, otherwise the owning
1070 // processor must have set the correct entry
1071 nonlocal_matrix->FillComplete(*column_space_map, matrix->RowMap());
1072 if (nonlocal_matrix_exporter.get() == nullptr)
1074 std::make_unique<Epetra_Export>(nonlocal_matrix->RowMap(),
1075 matrix->RowMap());
1076 ierr =
1078 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1079 ierr = matrix->FillComplete(*column_space_map, matrix->RowMap());
1080 nonlocal_matrix->PutScalar(0);
1081 }
1082 else
1083 ierr =
1084 matrix->GlobalAssemble(*column_space_map, matrix->RowMap(), true, mode);
1085
1086 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1087
1088 ierr = matrix->OptimizeStorage();
1089 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1090
1091 last_action = Zero;
1092
1093 compressed = true;
1094 }
1095
1096
1097
1098 void
1100 {
1101 // When we clear the matrix, reset
1102 // the pointer and generate an
1103 // empty matrix.
1105 std::make_unique<Epetra_Map>(0, 0, Utilities::Trilinos::comm_self());
1106 matrix = std::make_unique<Epetra_FECrsMatrix>(View, *column_space_map, 0);
1107 nonlocal_matrix.reset();
1109
1110 matrix->FillComplete();
1111
1112 compressed = true;
1113 }
1114
1115
1116
1117 void
1119 const TrilinosScalar new_diag_value)
1120 {
1121 Assert(matrix->Filled() == true, ExcMatrixNotCompressed());
1122
1123 // Only do this on the rows owned
1124 // locally on this processor.
1125 int local_row =
1126 matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
1127 if (local_row >= 0)
1128 {
1129 TrilinosScalar *values;
1130 int *col_indices;
1131 int num_entries;
1132 const int ierr =
1133 matrix->ExtractMyRowView(local_row, num_entries, values, col_indices);
1134 (void)ierr;
1135
1136 Assert(ierr == 0, ExcTrilinosError(ierr));
1137
1138 const std::ptrdiff_t diag_index =
1139 std::find(col_indices, col_indices + num_entries, local_row) -
1140 col_indices;
1141
1142 for (TrilinosWrappers::types::int_type j = 0; j < num_entries; ++j)
1143 if (diag_index != j || new_diag_value == 0)
1144 values[j] = 0.;
1145
1146 if (diag_index != num_entries)
1147 values[diag_index] = new_diag_value;
1148 }
1149 }
1150
1151
1152
1153 void
1155 const TrilinosScalar new_diag_value)
1156 {
1157 for (const auto row : rows)
1158 clear_row(row, new_diag_value);
1159 }
1160
1161
1162
1165 {
1166 // Extract local indices in
1167 // the matrix.
1168 int trilinos_i =
1169 matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
1170 trilinos_j =
1171 matrix->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
1172 TrilinosScalar value = 0.;
1173
1174 // If the data is not on the
1175 // present processor, we throw
1176 // an exception. This is one of
1177 // the two tiny differences to
1178 // the el(i,j) call, which does
1179 // not throw any assertions.
1180 if (trilinos_i == -1)
1181 {
1182 Assert(false,
1184 i, j, local_range().first, local_range().second - 1));
1185 }
1186 else
1187 {
1188 // Check whether the matrix has
1189 // already been transformed to local
1190 // indices.
1191 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1192
1193 // Prepare pointers for extraction
1194 // of a view of the row.
1195 int nnz_present = matrix->NumMyEntries(trilinos_i);
1196 int nnz_extracted;
1197 int *col_indices;
1198 TrilinosScalar *values;
1199
1200 // Generate the view and make
1201 // sure that we have not generated
1202 // an error.
1203 // TODO Check that col_indices are int and not long long
1204 int ierr = matrix->ExtractMyRowView(trilinos_i,
1205 nnz_extracted,
1206 values,
1207 col_indices);
1208 (void)ierr;
1209 Assert(ierr == 0, ExcTrilinosError(ierr));
1210
1211 Assert(nnz_present == nnz_extracted,
1212 ExcDimensionMismatch(nnz_present, nnz_extracted));
1213
1214 // Search the index where we
1215 // look for the value, and then
1216 // finally get it.
1217 const std::ptrdiff_t local_col_index =
1218 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1219 col_indices;
1220
1221 // This is actually the only
1222 // difference to the el(i,j)
1223 // function, which means that
1224 // we throw an exception in
1225 // this case instead of just
1226 // returning zero for an
1227 // element that is not present
1228 // in the sparsity pattern.
1229 if (local_col_index == nnz_present)
1230 {
1231 Assert(false, ExcInvalidIndex(i, j));
1232 }
1233 else
1234 value = values[local_col_index];
1235 }
1236
1237 return value;
1238 }
1239
1240
1241
1243 SparseMatrix::el(const size_type i, const size_type j) const
1244 {
1245 // Extract local indices in
1246 // the matrix.
1247 int trilinos_i =
1248 matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
1249 trilinos_j =
1250 matrix->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
1251 TrilinosScalar value = 0.;
1252
1253 // If the data is not on the
1254 // present processor, we can't
1255 // continue. Just print out zero
1256 // as discussed in the
1257 // documentation of this
1258 // function. if you want error
1259 // checking, use operator().
1260 if ((trilinos_i == -1) || (trilinos_j == -1))
1261 return 0.;
1262 else
1263 {
1264 // Check whether the matrix
1265 // already is transformed to
1266 // local indices.
1267 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1268
1269 // Prepare pointers for extraction
1270 // of a view of the row.
1271 int nnz_present = matrix->NumMyEntries(trilinos_i);
1272 int nnz_extracted;
1273 int *col_indices;
1274 TrilinosScalar *values;
1275
1276 // Generate the view and make
1277 // sure that we have not generated
1278 // an error.
1279 int ierr = matrix->ExtractMyRowView(trilinos_i,
1280 nnz_extracted,
1281 values,
1282 col_indices);
1283 (void)ierr;
1284 Assert(ierr == 0, ExcTrilinosError(ierr));
1285
1286 Assert(nnz_present == nnz_extracted,
1287 ExcDimensionMismatch(nnz_present, nnz_extracted));
1288
1289 // Search the index where we
1290 // look for the value, and then
1291 // finally get it.
1292 const std::ptrdiff_t local_col_index =
1293 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1294 col_indices;
1295
1296 // This is actually the only
1297 // difference to the () function
1298 // querying (i,j), where we throw an
1299 // exception instead of just
1300 // returning zero for an element
1301 // that is not present in the
1302 // sparsity pattern.
1303 if (local_col_index == nnz_present)
1304 value = 0;
1305 else
1306 value = values[local_col_index];
1307 }
1308
1309 return value;
1310 }
1311
1312
1313
1316 {
1317 Assert(m() == n(), ExcNotQuadratic());
1318
1319# ifdef DEBUG
1320 // use operator() in debug mode because
1321 // it checks if this is a valid element
1322 // (in parallel)
1323 return operator()(i, i);
1324# else
1325 // Trilinos doesn't seem to have a
1326 // more efficient way to access the
1327 // diagonal than by just using the
1328 // standard el(i,j) function.
1329 return el(i, i);
1330# endif
1331 }
1332
1333
1334
1335 unsigned int
1337 {
1338 Assert(row < m(), ExcInternalError());
1339
1340 // get a representation of the
1341 // present row
1342 int ncols = -1;
1343 int local_row =
1344 matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
1345 Assert((local_row >= 0), ExcAccessToNonlocalRow(row));
1346
1347 // on the processor who owns this
1348 // row, we'll have a non-negative
1349 // value.
1350 if (local_row >= 0)
1351 {
1352 int ierr = matrix->NumMyRowEntries(local_row, ncols);
1353 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1354 }
1355
1356 return static_cast<unsigned int>(ncols);
1357 }
1358
1359
1360
1361 void
1362 SparseMatrix::set(const std::vector<size_type> &row_indices,
1363 const std::vector<size_type> &col_indices,
1364 const FullMatrix<TrilinosScalar> &values,
1365 const bool elide_zero_values)
1366 {
1367 Assert(row_indices.size() == values.m(),
1368 ExcDimensionMismatch(row_indices.size(), values.m()));
1369 Assert(col_indices.size() == values.n(),
1370 ExcDimensionMismatch(col_indices.size(), values.n()));
1371
1372 for (size_type i = 0; i < row_indices.size(); ++i)
1373 set(row_indices[i],
1374 col_indices.size(),
1375 col_indices.data(),
1376 &values(i, 0),
1377 elide_zero_values);
1378 }
1379
1380
1381
1382 void
1384 const std::vector<size_type> &col_indices,
1385 const std::vector<TrilinosScalar> &values,
1386 const bool elide_zero_values)
1387 {
1388 Assert(col_indices.size() == values.size(),
1389 ExcDimensionMismatch(col_indices.size(), values.size()));
1390
1391 set(row,
1392 col_indices.size(),
1393 col_indices.data(),
1394 values.data(),
1395 elide_zero_values);
1396 }
1397
1398
1399
1400 template <>
1401 void
1402 SparseMatrix::set<TrilinosScalar>(const size_type row,
1403 const size_type n_cols,
1404 const size_type *col_indices,
1405 const TrilinosScalar *values,
1406 const bool elide_zero_values)
1407 {
1408 AssertIndexRange(row, this->m());
1409
1410 int ierr;
1411 if (last_action == Add)
1412 {
1413 ierr =
1414 matrix->GlobalAssemble(*column_space_map, matrix->RowMap(), true);
1415
1416 Assert(ierr == 0, ExcTrilinosError(ierr));
1417 }
1418
1419 last_action = Insert;
1420
1421 const TrilinosWrappers::types::int_type *col_index_ptr;
1422 const TrilinosScalar *col_value_ptr;
1423 const TrilinosWrappers::types::int_type trilinos_row = row;
1425
1426 boost::container::small_vector<TrilinosScalar, 200> local_value_array(
1427 elide_zero_values ? n_cols : 0);
1428 boost::container::small_vector<TrilinosWrappers::types::int_type, 200>
1429 local_index_array(elide_zero_values ? n_cols : 0);
1430
1431 // If we don't elide zeros, the pointers are already available... need to
1432 // cast to non-const pointers as that is the format taken by Trilinos (but
1433 // we will not modify const data)
1434 if (elide_zero_values == false)
1435 {
1436 col_index_ptr =
1437 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1438 col_indices);
1439 col_value_ptr = values;
1440 n_columns = n_cols;
1441 }
1442 else
1443 {
1444 // Otherwise, extract nonzero values in each row and get the
1445 // respective indices.
1446 col_index_ptr = local_index_array.data();
1447 col_value_ptr = local_value_array.data();
1448
1449 n_columns = 0;
1450 for (size_type j = 0; j < n_cols; ++j)
1451 {
1452 const double value = values[j];
1453 AssertIsFinite(value);
1454 if (value != 0)
1455 {
1456 local_index_array[n_columns] = col_indices[j];
1457 local_value_array[n_columns] = value;
1458 ++n_columns;
1459 }
1460 }
1461
1462 AssertIndexRange(n_columns, n_cols + 1);
1463 }
1464
1465
1466 // If the calling matrix owns the row to which we want to insert values,
1467 // we can directly call the Epetra_CrsMatrix input function, which is much
1468 // faster than the Epetra_FECrsMatrix function. We distinguish between two
1469 // cases: the first one is when the matrix is not filled (i.e., it is
1470 // possible to add new elements to the sparsity pattern), and the second
1471 // one is when the pattern is already fixed. In the former case, we add
1472 // the possibility to insert new values, and in the second we just replace
1473 // data.
1474 if (matrix->RowMap().MyGID(
1475 static_cast<TrilinosWrappers::types::int_type>(row)) == true)
1476 {
1477 if (matrix->Filled() == false)
1478 {
1479 ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(
1480 row, static_cast<int>(n_columns), col_value_ptr, col_index_ptr);
1481
1482 // When inserting elements, we do not want to create exceptions in
1483 // the case when inserting non-local data (since that's what we
1484 // want to do right now).
1485 if (ierr > 0)
1486 ierr = 0;
1487 }
1488 else
1489 ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row,
1490 n_columns,
1491 col_value_ptr,
1492 col_index_ptr);
1493 }
1494 else
1495 {
1496 // When we're at off-processor data, we have to stick with the
1497 // standard Insert/ReplaceGlobalValues function. Nevertheless, the way
1498 // we call it is the fastest one (any other will lead to repeated
1499 // allocation and deallocation of memory in order to call the function
1500 // we already use, which is very inefficient if writing one element at
1501 // a time).
1502 compressed = false;
1503
1504 if (matrix->Filled() == false)
1505 {
1506 ierr = matrix->InsertGlobalValues(1,
1507 &trilinos_row,
1508 n_columns,
1509 col_index_ptr,
1510 &col_value_ptr,
1511 Epetra_FECrsMatrix::ROW_MAJOR);
1512 if (ierr > 0)
1513 ierr = 0;
1514 }
1515 else
1516 ierr = matrix->ReplaceGlobalValues(1,
1517 &trilinos_row,
1518 n_columns,
1519 col_index_ptr,
1520 &col_value_ptr,
1521 Epetra_FECrsMatrix::ROW_MAJOR);
1522 // use the FECrsMatrix facilities for set even in the case when we
1523 // have explicitly set the off-processor rows because that only works
1524 // properly when adding elements, not when setting them (since we want
1525 // to only touch elements that have been set explicitly, and there is
1526 // no way on the receiving processor to identify them otherwise)
1527 }
1528
1529 Assert(ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
1530 AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
1531 }
1532
1533
1534
1535 void
1536 SparseMatrix::add(const std::vector<size_type> &indices,
1537 const FullMatrix<TrilinosScalar> &values,
1538 const bool elide_zero_values)
1539 {
1540 Assert(indices.size() == values.m(),
1541 ExcDimensionMismatch(indices.size(), values.m()));
1542 Assert(values.m() == values.n(), ExcNotQuadratic());
1543
1544 for (size_type i = 0; i < indices.size(); ++i)
1545 add(indices[i],
1546 indices.size(),
1547 indices.data(),
1548 &values(i, 0),
1549 elide_zero_values);
1550 }
1551
1552
1553
1554 void
1555 SparseMatrix::add(const std::vector<size_type> &row_indices,
1556 const std::vector<size_type> &col_indices,
1557 const FullMatrix<TrilinosScalar> &values,
1558 const bool elide_zero_values)
1559 {
1560 Assert(row_indices.size() == values.m(),
1561 ExcDimensionMismatch(row_indices.size(), values.m()));
1562 Assert(col_indices.size() == values.n(),
1563 ExcDimensionMismatch(col_indices.size(), values.n()));
1564
1565 for (size_type i = 0; i < row_indices.size(); ++i)
1566 add(row_indices[i],
1567 col_indices.size(),
1568 col_indices.data(),
1569 &values(i, 0),
1570 elide_zero_values);
1571 }
1572
1573
1574
1575 void
1577 const std::vector<size_type> &col_indices,
1578 const std::vector<TrilinosScalar> &values,
1579 const bool elide_zero_values)
1580 {
1581 Assert(col_indices.size() == values.size(),
1582 ExcDimensionMismatch(col_indices.size(), values.size()));
1583
1584 add(row,
1585 col_indices.size(),
1586 col_indices.data(),
1587 values.data(),
1588 elide_zero_values);
1589 }
1590
1591
1592
1593 void
1595 const size_type n_cols,
1596 const size_type *col_indices,
1597 const TrilinosScalar *values,
1598 const bool elide_zero_values,
1599 const bool /*col_indices_are_sorted*/)
1600 {
1601 AssertIndexRange(row, this->m());
1602 for (size_type n = 0; n < n_cols; ++n)
1603 AssertIndexRange(col_indices[n], this->n());
1604
1605 int ierr;
1606 if (last_action == Insert)
1607 {
1608 // TODO: this could lead to a dead lock when only one processor
1609 // calls GlobalAssemble.
1610 ierr =
1611 matrix->GlobalAssemble(*column_space_map, matrix->RowMap(), false);
1612
1613 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1614 }
1615
1616 last_action = Add;
1617
1618 const TrilinosWrappers::types::int_type *col_index_ptr;
1619 const TrilinosScalar *col_value_ptr;
1620 const TrilinosWrappers::types::int_type trilinos_row = row;
1622
1623 boost::container::small_vector<TrilinosScalar, 100> local_value_array(
1624 n_cols);
1625 boost::container::small_vector<TrilinosWrappers::types::int_type, 100>
1626 local_index_array(n_cols);
1627
1628 // If we don't elide zeros, the pointers are already available... need to
1629 // cast to non-const pointers as that is the format taken by Trilinos (but
1630 // we will not modify const data)
1631 if (elide_zero_values == false)
1632 {
1633 col_index_ptr =
1634 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1635 col_indices);
1636 col_value_ptr = values;
1637 n_columns = n_cols;
1638# ifdef DEBUG
1639 for (size_type j = 0; j < n_cols; ++j)
1640 AssertIsFinite(values[j]);
1641# endif
1642 }
1643 else
1644 {
1645 // Otherwise, extract nonzero values in each row and the corresponding
1646 // index.
1647 col_index_ptr = local_index_array.data();
1648 col_value_ptr = local_value_array.data();
1649
1650 n_columns = 0;
1651 for (size_type j = 0; j < n_cols; ++j)
1652 {
1653 const double value = values[j];
1654
1655 AssertIsFinite(value);
1656 if (value != 0)
1657 {
1658 local_index_array[n_columns] = col_indices[j];
1659 local_value_array[n_columns] = value;
1660 ++n_columns;
1661 }
1662 }
1663
1664 AssertIndexRange(n_columns, n_cols + 1);
1665 }
1666 // Exit early if there is nothing to do
1667 if (n_columns == 0)
1668 {
1669 return;
1670 }
1671
1672 // If the calling processor owns the row to which we want to add values, we
1673 // can directly call the Epetra_CrsMatrix input function, which is much
1674 // faster than the Epetra_FECrsMatrix function.
1675 if (matrix->RowMap().MyGID(
1676 static_cast<TrilinosWrappers::types::int_type>(row)) == true)
1677 {
1678 ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row,
1679 n_columns,
1680 col_value_ptr,
1681 col_index_ptr);
1682 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1683 }
1684 else if (nonlocal_matrix.get() != nullptr)
1685 {
1686 compressed = false;
1687 // this is the case when we have explicitly set the off-processor rows
1688 // and want to create a separate matrix object for them (to retain
1689 // thread-safety)
1690 Assert(nonlocal_matrix->RowMap().LID(
1691 static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1692 ExcMessage("Attempted to write into off-processor matrix row "
1693 "that has not be specified as being writable upon "
1694 "initialization"));
1695 ierr = nonlocal_matrix->SumIntoGlobalValues(row,
1696 n_columns,
1697 col_value_ptr,
1698 col_index_ptr);
1699 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1700 }
1701 else
1702 {
1703 // When we're at off-processor data, we have to stick with the
1704 // standard SumIntoGlobalValues function. Nevertheless, the way we
1705 // call it is the fastest one (any other will lead to repeated
1706 // allocation and deallocation of memory in order to call the function
1707 // we already use, which is very inefficient if writing one element at
1708 // a time).
1709 compressed = false;
1710
1711 ierr = matrix->SumIntoGlobalValues(1,
1712 &trilinos_row,
1713 n_columns,
1714 col_index_ptr,
1715 &col_value_ptr,
1716 Epetra_FECrsMatrix::ROW_MAJOR);
1717 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1718 }
1719
1720# ifdef DEBUG
1721 if (ierr > 0)
1722 {
1723 std::cout << "------------------------------------------" << std::endl;
1724 std::cout << "Got error " << ierr << " in row " << row << " of proc "
1725 << matrix->RowMap().Comm().MyPID()
1726 << " when trying to add the columns:" << std::endl;
1727 for (TrilinosWrappers::types::int_type i = 0; i < n_columns; ++i)
1728 std::cout << col_index_ptr[i] << " ";
1729 std::cout << std::endl << std::endl;
1730 std::cout << "Matrix row "
1731 << (matrix->RowMap().MyGID(
1732 static_cast<TrilinosWrappers::types::int_type>(row)) ==
1733 false ?
1734 "(nonlocal part)" :
1735 "")
1736 << " has the following indices:" << std::endl;
1737 std::vector<TrilinosWrappers::types::int_type> indices;
1738 const Epetra_CrsGraph *graph =
1739 (nonlocal_matrix.get() != nullptr &&
1740 matrix->RowMap().MyGID(
1741 static_cast<TrilinosWrappers::types::int_type>(row)) == false) ?
1742 &nonlocal_matrix->Graph() :
1743 &matrix->Graph();
1744
1745 indices.resize(graph->NumGlobalIndices(row));
1746 int n_indices = 0;
1747 graph->ExtractGlobalRowCopy(row,
1748 indices.size(),
1749 n_indices,
1750 indices.data());
1751 AssertDimension(n_indices, indices.size());
1752
1753 for (TrilinosWrappers::types::int_type i = 0; i < n_indices; ++i)
1754 std::cout << indices[i] << " ";
1755 std::cout << std::endl << std::endl;
1756 Assert(ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
1757 }
1758# endif
1759 Assert(ierr >= 0, ExcTrilinosError(ierr));
1760 }
1761
1762
1763
1764 SparseMatrix &
1766 {
1767 (void)d;
1769 compress(VectorOperation::unknown); // TODO: why do we do this? Should we
1770 // not check for is_compressed?
1771
1772 // As checked above, we are only allowed to use d==0.0, so pass
1773 // a constant zero (instead of a run-time value 'd' that *happens* to
1774 // have a zero value) to the underlying class in hopes that the compiler
1775 // can optimize this somehow.
1776 const int ierr = matrix->PutScalar(/*d=*/0.0);
1777 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1778
1779 if (nonlocal_matrix.get() != nullptr)
1780 {
1781 const int ierr = nonlocal_matrix->PutScalar(/*d=*/0.0);
1782 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1783 }
1784
1785 return *this;
1786 }
1787
1788
1789
1790 void
1792 {
1793 AssertDimension(rhs.m(), m());
1794 AssertDimension(rhs.n(), n());
1797 Assert(matrix->RowMap().SameAs(rhs.matrix->RowMap()),
1798 ExcMessage("Can only add matrices with same distribution of rows"));
1799 Assert(matrix->Filled() && rhs.matrix->Filled(),
1800 ExcMessage("Addition of matrices only allowed if matrices are "
1801 "filled, i.e., compress() has been called"));
1802
1803 const bool same_col_map = matrix->ColMap().SameAs(rhs.matrix->ColMap());
1804
1805 for (const auto row : locally_owned_range_indices())
1806 {
1807 const int row_local = matrix->RowMap().LID(
1808 static_cast<TrilinosWrappers::types::int_type>(row));
1809 Assert((row_local >= 0), ExcAccessToNonlocalRow(row));
1810
1811 // First get a view to the matrix columns of both matrices. Note that
1812 // the data is in local index spaces so we need to be careful not only
1813 // to compare column indices in case they are derived from the same
1814 // map.
1815 int n_entries, rhs_n_entries;
1816 TrilinosScalar *value_ptr, *rhs_value_ptr;
1817 int *index_ptr, *rhs_index_ptr;
1818 int ierr = rhs.matrix->ExtractMyRowView(row_local,
1819 rhs_n_entries,
1820 rhs_value_ptr,
1821 rhs_index_ptr);
1822 (void)ierr;
1823 Assert(ierr == 0, ExcTrilinosError(ierr));
1824
1825 ierr =
1826 matrix->ExtractMyRowView(row_local, n_entries, value_ptr, index_ptr);
1827 Assert(ierr == 0, ExcTrilinosError(ierr));
1828 bool expensive_checks = (n_entries != rhs_n_entries || !same_col_map);
1829 if (!expensive_checks)
1830 {
1831 // check if the column indices are the same. If yes, can simply
1832 // copy over the data.
1833 expensive_checks = std::memcmp(static_cast<void *>(index_ptr),
1834 static_cast<void *>(rhs_index_ptr),
1835 sizeof(int) * n_entries) != 0;
1836 if (!expensive_checks)
1837 for (int i = 0; i < n_entries; ++i)
1838 value_ptr[i] += rhs_value_ptr[i] * factor;
1839 }
1840 // Now to the expensive case where we need to check all column indices
1841 // against each other (transformed into global index space) and where
1842 // we need to make sure that all entries we are about to add into the
1843 // lhs matrix actually exist
1844 if (expensive_checks)
1845 {
1846 for (int i = 0; i < rhs_n_entries; ++i)
1847 {
1848 if (rhs_value_ptr[i] == 0.)
1849 continue;
1850 const TrilinosWrappers::types::int_type rhs_global_col =
1851 global_column_index(*rhs.matrix, rhs_index_ptr[i]);
1852 int local_col = matrix->ColMap().LID(rhs_global_col);
1853 int *local_index = Utilities::lower_bound(index_ptr,
1854 index_ptr + n_entries,
1855 local_col);
1856 Assert(local_index != index_ptr + n_entries &&
1857 *local_index == local_col,
1858 ExcMessage(
1859 "Adding the entries from the other matrix "
1860 "failed, because the sparsity pattern "
1861 "of that matrix includes more elements than the "
1862 "calling matrix, which is not allowed."));
1863 value_ptr[local_index - index_ptr] += factor * rhs_value_ptr[i];
1864 }
1865 }
1866 }
1867 }
1868
1869
1870
1871 void
1873 {
1874 // This only flips a flag that tells
1875 // Trilinos that any vmult operation
1876 // should be done with the
1877 // transpose. However, the matrix
1878 // structure is not reset.
1879 int ierr;
1880
1881 if (!matrix->UseTranspose())
1882 {
1883 ierr = matrix->SetUseTranspose(true);
1884 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1885 }
1886 else
1887 {
1888 ierr = matrix->SetUseTranspose(false);
1889 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1890 }
1891 }
1892
1893
1894
1895 SparseMatrix &
1897 {
1898 const int ierr = matrix->Scale(a);
1899 Assert(ierr == 0, ExcTrilinosError(ierr));
1900 (void)ierr; // removes -Wunused-variable in optimized mode
1901
1902 return *this;
1903 }
1904
1905
1906
1907 SparseMatrix &
1909 {
1910 Assert(a != 0, ExcDivideByZero());
1911
1912 const TrilinosScalar factor = 1. / a;
1913
1914 const int ierr = matrix->Scale(factor);
1915 Assert(ierr == 0, ExcTrilinosError(ierr));
1916 (void)ierr; // removes -Wunused-variable in optimized mode
1917
1918 return *this;
1919 }
1920
1921
1922
1925 {
1926 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1927 return matrix->NormOne();
1928 }
1929
1930
1931
1934 {
1935 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1936 return matrix->NormInf();
1937 }
1938
1939
1940
1943 {
1944 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1945 return matrix->NormFrobenius();
1946 }
1947
1948
1949
1950 namespace internal
1951 {
1952 namespace SparseMatrixImplementation
1953 {
1954 template <typename VectorType>
1955 inline void
1956 check_vector_map_equality(const Epetra_CrsMatrix &,
1957 const VectorType &,
1958 const VectorType &)
1959 {}
1960
1961 inline void
1962 check_vector_map_equality(const Epetra_CrsMatrix &m,
1965 {
1966 Assert(in.trilinos_partitioner().SameAs(m.DomainMap()) == true,
1967 ExcMessage("The column partitioning of a matrix does not match "
1968 "the partitioning of a vector you are trying to "
1969 "multiply it with. Are you multiplying the "
1970 "matrix with a vector that has ghost elements?"));
1971 Assert(out.trilinos_partitioner().SameAs(m.RangeMap()) == true,
1972 ExcMessage("The row partitioning of a matrix does not match "
1973 "the partitioning of a vector you are trying to "
1974 "put the result of a matrix-vector product in. "
1975 "Are you trying to put the product of the "
1976 "matrix with a vector into a vector that has "
1977 "ghost elements?"));
1978 (void)m;
1979 (void)in;
1980 (void)out;
1981 }
1982 } // namespace SparseMatrixImplementation
1983 } // namespace internal
1984
1985
1986 template <typename VectorType>
1987 std::enable_if_t<
1988 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1989 SparseMatrix::vmult(VectorType &dst, const VectorType &src) const
1990 {
1991 Assert(&src != &dst, ExcSourceEqualsDestination());
1992 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1993 (void)src;
1994 (void)dst;
1995
1997 src,
1998 dst);
1999 const size_type dst_local_size = internal::end(dst) - internal::begin(dst);
2000 AssertDimension(dst_local_size, matrix->RangeMap().NumMyPoints());
2001 const size_type src_local_size = internal::end(src) - internal::begin(src);
2002 AssertDimension(src_local_size, matrix->DomainMap().NumMyPoints());
2003
2004 Epetra_MultiVector tril_dst(
2005 View, matrix->RangeMap(), internal::begin(dst), dst_local_size, 1);
2006 Epetra_MultiVector tril_src(View,
2007 matrix->DomainMap(),
2008 const_cast<TrilinosScalar *>(
2009 internal::begin(src)),
2010 src_local_size,
2011 1);
2012
2013 const int ierr = matrix->Multiply(false, tril_src, tril_dst);
2014 Assert(ierr == 0, ExcTrilinosError(ierr));
2015 (void)ierr; // removes -Wunused-variable in optimized mode
2016 }
2017
2018
2019
2020 template <typename VectorType>
2021 std::enable_if_t<
2022 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2023 SparseMatrix::vmult(VectorType & /*dst*/, const VectorType & /*src*/) const
2024 {
2026 }
2027
2028
2029
2030 template <typename VectorType>
2031 std::enable_if_t<
2032 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2033 SparseMatrix::Tvmult(VectorType &dst, const VectorType &src) const
2034 {
2035 Assert(&src != &dst, ExcSourceEqualsDestination());
2036 Assert(matrix->Filled(), ExcMatrixNotCompressed());
2037
2039 dst,
2040 src);
2041 const size_type dst_local_size = internal::end(dst) - internal::begin(dst);
2042 AssertDimension(dst_local_size, matrix->DomainMap().NumMyPoints());
2043 const size_type src_local_size = internal::end(src) - internal::begin(src);
2044 AssertDimension(src_local_size, matrix->RangeMap().NumMyPoints());
2045
2046 Epetra_MultiVector tril_dst(
2047 View, matrix->DomainMap(), internal::begin(dst), dst_local_size, 1);
2048 Epetra_MultiVector tril_src(View,
2049 matrix->RangeMap(),
2050 const_cast<double *>(internal::begin(src)),
2051 src_local_size,
2052 1);
2053
2054 const int ierr = matrix->Multiply(true, tril_src, tril_dst);
2055 Assert(ierr == 0, ExcTrilinosError(ierr));
2056 (void)ierr; // removes -Wunused-variable in optimized mode
2057 }
2058
2059
2060
2061 template <typename VectorType>
2062 std::enable_if_t<
2063 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2064 SparseMatrix::Tvmult(VectorType & /*dst*/, const VectorType & /*src*/) const
2065 {
2067 }
2068
2069
2070
2071 template <typename VectorType>
2072 void
2073 SparseMatrix::vmult_add(VectorType &dst, const VectorType &src) const
2074 {
2075 Assert(&src != &dst, ExcSourceEqualsDestination());
2076
2077 // Reinit a temporary vector with fast argument set, which does not
2078 // overwrite the content (to save time).
2079 VectorType tmp_vector;
2080 tmp_vector.reinit(dst, true);
2081 vmult(tmp_vector, src);
2082 dst += tmp_vector;
2083 }
2084
2085
2086
2087 template <typename VectorType>
2088 void
2089 SparseMatrix::Tvmult_add(VectorType &dst, const VectorType &src) const
2090 {
2091 Assert(&src != &dst, ExcSourceEqualsDestination());
2092
2093 // Reinit a temporary vector with fast argument set, which does not
2094 // overwrite the content (to save time).
2095 VectorType tmp_vector;
2096 tmp_vector.reinit(dst, true);
2097 Tvmult(tmp_vector, src);
2098 dst += tmp_vector;
2099 }
2100
2101
2102
2105 {
2106 AssertDimension(m(), v.size());
2107 Assert(matrix->RowMap().SameAs(matrix->DomainMap()), ExcNotQuadratic());
2108
2109 MPI::Vector temp_vector;
2110 temp_vector.reinit(v, true);
2111
2112 vmult(temp_vector, v);
2113 return temp_vector * v;
2114 }
2115
2116
2117
2120 const MPI::Vector &v) const
2121 {
2122 AssertDimension(m(), u.size());
2123 AssertDimension(m(), v.size());
2124 Assert(matrix->RowMap().SameAs(matrix->DomainMap()), ExcNotQuadratic());
2125
2126 MPI::Vector temp_vector;
2127 temp_vector.reinit(v, true);
2128
2129 vmult(temp_vector, v);
2130 return u * temp_vector;
2131 }
2132
2133
2134
2135 namespace internals
2136 {
2137 void
2138 perform_mmult(const SparseMatrix &inputleft,
2139 const SparseMatrix &inputright,
2140 SparseMatrix &result,
2141 const MPI::Vector &V,
2142 const bool transpose_left)
2143 {
2144 const bool use_vector = (V.size() == inputright.m() ? true : false);
2145 if (transpose_left == false)
2146 {
2147 Assert(inputleft.n() == inputright.m(),
2148 ExcDimensionMismatch(inputleft.n(), inputright.m()));
2149 Assert(inputleft.trilinos_matrix().DomainMap().SameAs(
2150 inputright.trilinos_matrix().RangeMap()),
2151 ExcMessage("Parallel partitioning of A and B does not fit."));
2152 }
2153 else
2154 {
2155 Assert(inputleft.m() == inputright.m(),
2156 ExcDimensionMismatch(inputleft.m(), inputright.m()));
2157 Assert(inputleft.trilinos_matrix().RangeMap().SameAs(
2158 inputright.trilinos_matrix().RangeMap()),
2159 ExcMessage("Parallel partitioning of A and B does not fit."));
2160 }
2161
2162 result.clear();
2163
2164 // create a suitable operator B: in case
2165 // we do not use a vector, all we need to
2166 // do is to set the pointer. Otherwise,
2167 // we insert the data from B, but
2168 // multiply each row with the respective
2169 // vector element.
2170 Teuchos::RCP<Epetra_CrsMatrix> mod_B;
2171 if (use_vector == false)
2172 {
2173 mod_B = Teuchos::rcp(const_cast<Epetra_CrsMatrix *>(
2174 &inputright.trilinos_matrix()),
2175 false);
2176 }
2177 else
2178 {
2179 mod_B = Teuchos::rcp(
2180 new Epetra_CrsMatrix(Copy, inputright.trilinos_sparsity_pattern()),
2181 true);
2182 mod_B->FillComplete(inputright.trilinos_matrix().DomainMap(),
2183 inputright.trilinos_matrix().RangeMap());
2184 Assert(inputright.local_range() == V.local_range(),
2185 ExcMessage("Parallel distribution of matrix B and vector V "
2186 "does not match."));
2187
2188 const int local_N = inputright.local_size();
2189 for (int i = 0; i < local_N; ++i)
2190 {
2191 int N_entries = -1;
2192 double *new_data, *B_data;
2193 mod_B->ExtractMyRowView(i, N_entries, new_data);
2194 inputright.trilinos_matrix().ExtractMyRowView(i,
2195 N_entries,
2196 B_data);
2197 double value = V.trilinos_vector()[0][i];
2198 for (TrilinosWrappers::types::int_type j = 0; j < N_entries; ++j)
2199 new_data[j] = value * B_data[j];
2200 }
2201 }
2202
2203
2204 SparseMatrix tmp_result(transpose_left ?
2205 inputleft.locally_owned_domain_indices() :
2206 inputleft.locally_owned_range_indices(),
2207 inputright.locally_owned_domain_indices(),
2208 inputleft.get_mpi_communicator());
2209
2210# ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT
2211 EpetraExt::MatrixMatrix::Multiply(inputleft.trilinos_matrix(),
2212 transpose_left,
2213 *mod_B,
2214 false,
2215 const_cast<Epetra_CrsMatrix &>(
2216 tmp_result.trilinos_matrix()));
2217# else
2218 Assert(false,
2219 ExcMessage("This function requires that the Trilinos "
2220 "installation found while running the deal.II "
2221 "CMake scripts contains the optional Trilinos "
2222 "package 'EpetraExt'. However, this optional "
2223 "part of Trilinos was not found."));
2224# endif
2225 result.reinit(tmp_result.trilinos_matrix());
2226 }
2227 } // namespace internals
2228
2229
2230 void
2232 const SparseMatrix &B,
2233 const MPI::Vector &V) const
2234 {
2235 internals::perform_mmult(*this, B, C, V, false);
2236 }
2237
2238
2239
2240 void
2242 const SparseMatrix &B,
2243 const MPI::Vector &V) const
2244 {
2245 internals::perform_mmult(*this, B, C, V, true);
2246 }
2247
2248
2249
2250 void
2255
2256
2257
2258 // As of now, no particularly neat
2259 // output is generated in case of
2260 // multiple processors.
2261 void
2262 SparseMatrix::print(std::ostream &out,
2263 const bool print_detailed_trilinos_information) const
2264 {
2265 if (print_detailed_trilinos_information == true)
2266 out << *matrix;
2267 else
2268 {
2269 double *values;
2270 int *indices;
2271 int num_entries;
2272
2273 for (int i = 0; i < matrix->NumMyRows(); ++i)
2274 {
2275 const int ierr =
2276 matrix->ExtractMyRowView(i, num_entries, values, indices);
2277 (void)ierr;
2278 Assert(ierr == 0, ExcTrilinosError(ierr));
2279
2280 for (TrilinosWrappers::types::int_type j = 0; j < num_entries; ++j)
2282 << ","
2284 << ") " << values[j] << std::endl;
2285 }
2286 }
2287
2288 AssertThrow(out.fail() == false, ExcIO());
2289 }
2290
2291
2292
2295 {
2296 size_type static_memory =
2297 sizeof(*this) + sizeof(*matrix) + sizeof(*matrix->Graph().DataPtr());
2298 return (
2300 matrix->NumMyNonzeros() +
2301 sizeof(int) * local_size() + static_memory);
2302 }
2303
2304
2305
2306 MPI_Comm
2308 {
2309 const Epetra_MpiComm *mpi_comm =
2310 dynamic_cast<const Epetra_MpiComm *>(&matrix->RangeMap().Comm());
2311 Assert(mpi_comm != nullptr, ExcInternalError());
2312 return mpi_comm->Comm();
2313 }
2314} // namespace TrilinosWrappers
2315
2316
2317namespace TrilinosWrappers
2318{
2319 namespace internal
2320 {
2321 namespace LinearOperatorImplementation
2322 {
2324
2326 : use_transpose(false)
2327 , communicator(MPI_COMM_SELF)
2328 , domain_map(IndexSet().make_trilinos_map(communicator.Comm()))
2329 , range_map(IndexSet().make_trilinos_map(communicator.Comm()))
2330 {
2331 vmult = [](Range &, const Domain &) {
2332 Assert(false,
2333 ExcMessage("Uninitialized TrilinosPayload::vmult called "
2334 "(Default constructor)"));
2335 };
2336
2337 Tvmult = [](Domain &, const Range &) {
2338 Assert(false,
2339 ExcMessage("Uninitialized TrilinosPayload::Tvmult called "
2340 "(Default constructor)"));
2341 };
2342
2343 inv_vmult = [](Domain &, const Range &) {
2344 Assert(false,
2345 ExcMessage("Uninitialized TrilinosPayload::inv_vmult called "
2346 "(Default constructor)"));
2347 };
2348
2349 inv_Tvmult = [](Range &, const Domain &) {
2350 Assert(false,
2351 ExcMessage("Uninitialized TrilinosPayload::inv_Tvmult called "
2352 "(Default constructor)"));
2353 };
2354 }
2355
2356
2357
2359 const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2360 const TrilinosWrappers::SparseMatrix &matrix)
2361 : TrilinosPayload(const_cast<Epetra_CrsMatrix &>(
2362 matrix.trilinos_matrix()),
2363 /*op_supports_inverse_operations = */ false,
2364 matrix_exemplar.trilinos_matrix().UseTranspose(),
2365 matrix_exemplar.get_mpi_communicator(),
2366 matrix_exemplar.locally_owned_domain_indices(),
2367 matrix_exemplar.locally_owned_range_indices())
2368 {}
2369
2370
2371
2373 const TrilinosPayload &payload_exemplar,
2374 const TrilinosWrappers::SparseMatrix &matrix)
2375
2376 : TrilinosPayload(const_cast<Epetra_CrsMatrix &>(
2377 matrix.trilinos_matrix()),
2378 /*op_supports_inverse_operations = */ false,
2379 payload_exemplar.UseTranspose(),
2380 payload_exemplar.get_mpi_communicator(),
2381 payload_exemplar.locally_owned_domain_indices(),
2382 payload_exemplar.locally_owned_range_indices())
2383 {}
2384
2385
2386
2388 const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2389 const TrilinosWrappers::PreconditionBase &preconditioner)
2390 : TrilinosPayload(preconditioner.trilinos_operator(),
2391 /*op_supports_inverse_operations = */ true,
2392 matrix_exemplar.trilinos_matrix().UseTranspose(),
2393 matrix_exemplar.get_mpi_communicator(),
2394 matrix_exemplar.locally_owned_domain_indices(),
2395 matrix_exemplar.locally_owned_range_indices())
2396 {}
2397
2398
2399
2401 const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2402 const TrilinosWrappers::PreconditionBase &preconditioner)
2404 preconditioner.trilinos_operator(),
2405 /*op_supports_inverse_operations = */ true,
2406 preconditioner_exemplar.trilinos_operator().UseTranspose(),
2407 preconditioner_exemplar.get_mpi_communicator(),
2408 preconditioner_exemplar.locally_owned_domain_indices(),
2409 preconditioner_exemplar.locally_owned_range_indices())
2410 {}
2411
2412
2413
2415 const TrilinosPayload &payload_exemplar,
2416 const TrilinosWrappers::PreconditionBase &preconditioner)
2417 : TrilinosPayload(preconditioner.trilinos_operator(),
2418 /*op_supports_inverse_operations = */ true,
2419 payload_exemplar.UseTranspose(),
2420 payload_exemplar.get_mpi_communicator(),
2421 payload_exemplar.locally_owned_domain_indices(),
2422 payload_exemplar.locally_owned_range_indices())
2423 {}
2424
2425
2426
2428 : vmult(payload.vmult)
2429 , Tvmult(payload.Tvmult)
2430 , inv_vmult(payload.inv_vmult)
2431 , inv_Tvmult(payload.inv_Tvmult)
2432 , use_transpose(payload.use_transpose)
2433 , communicator(payload.communicator)
2434 , domain_map(payload.domain_map)
2435 , range_map(payload.range_map)
2436 {}
2437
2438
2439
2440 // Composite copy constructor
2441 // This is required for PackagedOperations
2443 const TrilinosPayload &second_op)
2444 : use_transpose(false)
2445 , // The combination of operators provides the exact
2446 // definition of the operation
2447 communicator(first_op.communicator)
2448 , domain_map(second_op.domain_map)
2449 , range_map(first_op.range_map)
2450 {}
2451
2452
2453
2456 {
2457 TrilinosPayload return_op(*this);
2458
2459 return_op.vmult = [](Range &tril_dst, const Range &tril_src) {
2460 tril_dst = tril_src;
2461 };
2462
2463 return_op.Tvmult = [](Range &tril_dst, const Range &tril_src) {
2464 tril_dst = tril_src;
2465 };
2466
2467 return_op.inv_vmult = [](Range &tril_dst, const Range &tril_src) {
2468 tril_dst = tril_src;
2469 };
2470
2471 return_op.inv_Tvmult = [](Range &tril_dst, const Range &tril_src) {
2472 tril_dst = tril_src;
2473 };
2474
2475 return return_op;
2476 }
2477
2478
2479
2482 {
2483 TrilinosPayload return_op(*this);
2484
2485 return_op.vmult = [](Range &tril_dst, const Domain &) {
2486 const int ierr = tril_dst.PutScalar(0.0);
2487
2488 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2489 };
2490
2491 return_op.Tvmult = [](Domain &tril_dst, const Range &) {
2492 const int ierr = tril_dst.PutScalar(0.0);
2493
2494 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2495 };
2496
2497 return_op.inv_vmult = [](Domain &tril_dst, const Range &) {
2498 AssertThrow(false,
2499 ExcMessage("Cannot compute inverse of null operator"));
2500
2501 const int ierr = tril_dst.PutScalar(0.0);
2502
2503 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2504 };
2505
2506 return_op.inv_Tvmult = [](Range &tril_dst, const Domain &) {
2507 AssertThrow(false,
2508 ExcMessage("Cannot compute inverse of null operator"));
2509
2510 const int ierr = tril_dst.PutScalar(0.0);
2511
2512 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2513 };
2514
2515 return return_op;
2516 }
2517
2518
2519
2522 {
2523 TrilinosPayload return_op(*this);
2524 return_op.transpose();
2525 return return_op;
2526 }
2527
2528
2529
2530 IndexSet
2535
2536
2537
2538 IndexSet
2543
2544
2545
2546 MPI_Comm
2548 {
2549 return communicator.Comm();
2550 }
2551
2552
2553
2554 void
2559
2560
2561
2562 bool
2564 {
2565 return use_transpose;
2566 }
2567
2568
2569
2570 int
2572 {
2574 {
2576 std::swap(domain_map, range_map);
2577 std::swap(vmult, Tvmult);
2578 std::swap(inv_vmult, inv_Tvmult);
2579 }
2580 return 0;
2581 }
2582
2583
2584
2585 int
2587 {
2588 // The transposedness of the operations is taken care of
2589 // when we hit the transpose flag.
2590 vmult(Y, X);
2591 return 0;
2592 }
2593
2594
2595
2596 int
2598 {
2599 // The transposedness of the operations is taken care of
2600 // when we hit the transpose flag.
2601 inv_vmult(X, Y);
2602 return 0;
2603 }
2604
2605
2606
2607 const char *
2609 {
2610 return "TrilinosPayload";
2611 }
2612
2613
2614
2615 const Epetra_Comm &
2617 {
2618 return communicator;
2619 }
2620
2621
2622
2623 const Epetra_Map &
2625 {
2626 return domain_map;
2627 }
2628
2629
2630
2631 const Epetra_Map &
2633 {
2634 return range_map;
2635 }
2636
2637
2638
2639 bool
2641 {
2642 return false;
2643 }
2644
2645
2646
2647 double
2649 {
2651 return 0.0;
2652 }
2653
2654
2655
2658 const TrilinosPayload &second_op)
2659 {
2660 using Domain = typename TrilinosPayload::Domain;
2661 using Range = typename TrilinosPayload::Range;
2662 using Intermediate = typename TrilinosPayload::VectorType;
2663 using GVMVectorType = TrilinosWrappers::MPI::Vector;
2664
2666 second_op.locally_owned_domain_indices(),
2667 ExcMessage(
2668 "Operators are set to work on incompatible IndexSets."));
2670 second_op.locally_owned_range_indices(),
2671 ExcMessage(
2672 "Operators are set to work on incompatible IndexSets."));
2673
2674 TrilinosPayload return_op(first_op, second_op);
2675
2676 // Capture by copy so the payloads are always valid
2677 return_op.vmult = [first_op, second_op](Range &tril_dst,
2678 const Domain &tril_src) {
2679 // Duplicated from LinearOperator::operator*
2680 // TODO: Template the constructor on GrowingVectorMemory vector type?
2682 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2683
2684 // Initialize intermediate vector
2685 const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2686 i->reinit(IndexSet(first_op_init_map),
2687 first_op.get_mpi_communicator(),
2688 /*bool omit_zeroing_entries =*/true);
2689
2690 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2691 const size_type i_local_size = i->end() - i->begin();
2692 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2693 const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2694 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2695 (void)second_op_init_map;
2696 Intermediate tril_int(View,
2697 first_op_init_map,
2698 const_cast<TrilinosScalar *>(i->begin()),
2699 i_local_size,
2700 1);
2701
2702 // These operators may themselves be transposed or not, so we let them
2703 // decide what the intended outcome is
2704 second_op.Apply(tril_src, tril_int);
2705 first_op.Apply(tril_src, tril_dst);
2706 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2707 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2708 };
2709
2710 return_op.Tvmult = [first_op, second_op](Domain &tril_dst,
2711 const Range &tril_src) {
2712 // Duplicated from LinearOperator::operator*
2713 // TODO: Template the constructor on GrowingVectorMemory vector type?
2715 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2716
2717 // These operators may themselves be transposed or not, so we let them
2718 // decide what the intended outcome is
2719 // We must first transpose the operators to get the right IndexSets
2720 // for the input, intermediate and result vectors
2721 const_cast<TrilinosPayload &>(first_op).transpose();
2722 const_cast<TrilinosPayload &>(second_op).transpose();
2723
2724 // Initialize intermediate vector
2725 const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2726 i->reinit(IndexSet(first_op_init_map),
2727 first_op.get_mpi_communicator(),
2728 /*bool omit_zeroing_entries =*/true);
2729
2730 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2731 const size_type i_local_size = i->end() - i->begin();
2732 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2733 const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2734 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2735 (void)second_op_init_map;
2736 Intermediate tril_int(View,
2737 first_op_init_map,
2738 const_cast<TrilinosScalar *>(i->begin()),
2739 i_local_size,
2740 1);
2741
2742 // These operators may themselves be transposed or not, so we let them
2743 // decide what the intended outcome is
2744 second_op.Apply(tril_src, tril_int);
2745 first_op.Apply(tril_src, tril_dst);
2746 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2747 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2748
2749 // Reset transpose flag
2750 const_cast<TrilinosPayload &>(first_op).transpose();
2751 const_cast<TrilinosPayload &>(second_op).transpose();
2752 };
2753
2754 return_op.inv_vmult = [first_op, second_op](Domain &tril_dst,
2755 const Range &tril_src) {
2756 // Duplicated from LinearOperator::operator*
2757 // TODO: Template the constructor on GrowingVectorMemory vector type?
2759 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2760
2761 // Initialize intermediate vector
2762 const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2763 i->reinit(IndexSet(first_op_init_map),
2764 first_op.get_mpi_communicator(),
2765 /*bool omit_zeroing_entries =*/true);
2766
2767 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2768 const size_type i_local_size = i->end() - i->begin();
2769 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2770 const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2771 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2772 (void)second_op_init_map;
2773 Intermediate tril_int(View,
2774 first_op_init_map,
2775 const_cast<TrilinosScalar *>(i->begin()),
2776 i_local_size,
2777 1);
2778
2779 // These operators may themselves be transposed or not, so we let them
2780 // decide what the intended outcome is
2781 second_op.ApplyInverse(tril_src, tril_int);
2782 first_op.ApplyInverse(tril_src, tril_dst);
2783 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2784 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2785 };
2786
2787 return_op.inv_Tvmult = [first_op, second_op](Range &tril_dst,
2788 const Domain &tril_src) {
2789 // Duplicated from LinearOperator::operator*
2790 // TODO: Template the constructor on GrowingVectorMemory vector type?
2792 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2793
2794 // These operators may themselves be transposed or not, so we let them
2795 // decide what the intended outcome is
2796 // We must first transpose the operators to get the right IndexSets
2797 // for the input, intermediate and result vectors
2798 const_cast<TrilinosPayload &>(first_op).transpose();
2799 const_cast<TrilinosPayload &>(second_op).transpose();
2800
2801 // Initialize intermediate vector
2802 const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2803 i->reinit(IndexSet(first_op_init_map),
2804 first_op.get_mpi_communicator(),
2805 /*bool omit_zeroing_entries =*/true);
2806
2807 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2808 const size_type i_local_size = i->end() - i->begin();
2809 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2810 const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2811 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2812 (void)second_op_init_map;
2813 Intermediate tril_int(View,
2814 first_op_init_map,
2815 const_cast<TrilinosScalar *>(i->begin()),
2816 i_local_size,
2817 1);
2818
2819 // These operators may themselves be transposed or not, so we let them
2820 // decide what the intended outcome is
2821 second_op.ApplyInverse(tril_src, tril_int);
2822 first_op.ApplyInverse(tril_src, tril_dst);
2823 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2824 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2825
2826 // Reset transpose flag
2827 const_cast<TrilinosPayload &>(first_op).transpose();
2828 const_cast<TrilinosPayload &>(second_op).transpose();
2829 };
2830
2831 return return_op;
2832 }
2833
2834
2835
2836 TrilinosPayload
2838 const TrilinosPayload &second_op)
2839 {
2840 using Domain = typename TrilinosPayload::Domain;
2841 using Range = typename TrilinosPayload::Range;
2842 using Intermediate = typename TrilinosPayload::VectorType;
2843 using GVMVectorType = TrilinosWrappers::MPI::Vector;
2844
2846 second_op.locally_owned_range_indices(),
2847 ExcMessage(
2848 "Operators are set to work on incompatible IndexSets."));
2849
2850 TrilinosPayload return_op(first_op, second_op);
2851
2852 // Capture by copy so the payloads are always valid
2853 return_op.vmult = [first_op, second_op](Range &tril_dst,
2854 const Domain &tril_src) {
2855 // Duplicated from LinearOperator::operator*
2856 // TODO: Template the constructor on GrowingVectorMemory vector type?
2858 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2859
2860 // Initialize intermediate vector
2861 const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2862 i->reinit(IndexSet(first_op_init_map),
2863 first_op.get_mpi_communicator(),
2864 /*bool omit_zeroing_entries =*/true);
2865
2866 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2867 const size_type i_local_size = i->end() - i->begin();
2868 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2869 const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2870 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2871 (void)second_op_init_map;
2872 Intermediate tril_int(View,
2873 first_op_init_map,
2874 const_cast<TrilinosScalar *>(i->begin()),
2875 i_local_size,
2876 1);
2877
2878 // These operators may themselves be transposed or not, so we let them
2879 // decide what the intended outcome is
2880 second_op.Apply(tril_src, tril_int);
2881 first_op.Apply(tril_int, tril_dst);
2882 };
2883
2884 return_op.Tvmult = [first_op, second_op](Domain &tril_dst,
2885 const Range &tril_src) {
2886 // Duplicated from LinearOperator::operator*
2887 // TODO: Template the constructor on GrowingVectorMemory vector type?
2889 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2890
2891 // These operators may themselves be transposed or not, so we let them
2892 // decide what the intended outcome is
2893 // We must first transpose the operators to get the right IndexSets
2894 // for the input, intermediate and result vectors
2895 const_cast<TrilinosPayload &>(first_op).transpose();
2896 const_cast<TrilinosPayload &>(second_op).transpose();
2897
2898 // Initialize intermediate vector
2899 const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2900 i->reinit(IndexSet(first_op_init_map),
2901 first_op.get_mpi_communicator(),
2902 /*bool omit_zeroing_entries =*/true);
2903
2904 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2905 const size_type i_local_size = i->end() - i->begin();
2906 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2907 const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2908 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2909 (void)second_op_init_map;
2910 Intermediate tril_int(View,
2911 first_op_init_map,
2912 const_cast<TrilinosScalar *>(i->begin()),
2913 i_local_size,
2914 1);
2915
2916 // Apply the operators in the reverse order to vmult
2917 first_op.Apply(tril_src, tril_int);
2918 second_op.Apply(tril_int, tril_dst);
2919
2920 // Reset transpose flag
2921 const_cast<TrilinosPayload &>(first_op).transpose();
2922 const_cast<TrilinosPayload &>(second_op).transpose();
2923 };
2924
2925 return_op.inv_vmult = [first_op, second_op](Domain &tril_dst,
2926 const Range &tril_src) {
2927 // Duplicated from LinearOperator::operator*
2928 // TODO: Template the constructor on GrowingVectorMemory vector type?
2930 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2931
2932 // Initialize intermediate vector
2933 const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2934 i->reinit(IndexSet(first_op_init_map),
2935 first_op.get_mpi_communicator(),
2936 /*bool omit_zeroing_entries =*/true);
2937
2938 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2939 const size_type i_local_size = i->end() - i->begin();
2940 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2941 const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2942 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2943 (void)second_op_init_map;
2944 Intermediate tril_int(View,
2945 first_op_init_map,
2946 const_cast<TrilinosScalar *>(i->begin()),
2947 i_local_size,
2948 1);
2949
2950 // Apply the operators in the reverse order to vmult
2951 // and the same order as Tvmult
2952 first_op.ApplyInverse(tril_src, tril_int);
2953 second_op.ApplyInverse(tril_int, tril_dst);
2954 };
2955
2956 return_op.inv_Tvmult = [first_op, second_op](Range &tril_dst,
2957 const Domain &tril_src) {
2958 // Duplicated from LinearOperator::operator*
2959 // TODO: Template the constructor on GrowingVectorMemory vector type?
2961 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2962
2963 // These operators may themselves be transposed or not, so we let them
2964 // decide what the intended outcome is
2965 // We must first transpose the operators to get the right IndexSets
2966 // for the input, intermediate and result vectors
2967 const_cast<TrilinosPayload &>(first_op).transpose();
2968 const_cast<TrilinosPayload &>(second_op).transpose();
2969
2970 // Initialize intermediate vector
2971 const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2972 i->reinit(IndexSet(first_op_init_map),
2973 first_op.get_mpi_communicator(),
2974 /*bool omit_zeroing_entries =*/true);
2975
2976 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2977 const size_type i_local_size = i->end() - i->begin();
2978 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2979 const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2980 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2981 (void)second_op_init_map;
2982 Intermediate tril_int(View,
2983 first_op_init_map,
2984 const_cast<TrilinosScalar *>(i->begin()),
2985 i_local_size,
2986 1);
2987
2988 // These operators may themselves be transposed or not, so we let them
2989 // decide what the intended outcome is
2990 // Apply the operators in the reverse order to Tvmult
2991 // and the same order as vmult
2992 second_op.ApplyInverse(tril_src, tril_int);
2993 first_op.ApplyInverse(tril_int, tril_dst);
2994
2995 // Reset transpose flag
2996 const_cast<TrilinosPayload &>(first_op).transpose();
2997 const_cast<TrilinosPayload &>(second_op).transpose();
2998 };
2999
3000 return return_op;
3001 }
3002
3003 } // namespace LinearOperatorImplementation
3004 } /* namespace internal */
3005} /* namespace TrilinosWrappers */
3006
3007
3008
3009// explicit instantiations
3010# include "trilinos_sparse_matrix.inst"
3011
3012# ifndef DOXYGEN
3013// TODO: put these instantiations into generic file
3014namespace TrilinosWrappers
3015{
3016 template void
3017 SparseMatrix::reinit(const ::SparsityPattern &);
3018
3019 template void
3021
3022 template void
3024 const IndexSet &,
3025 const ::SparsityPattern &,
3026 const MPI_Comm,
3027 const bool);
3028
3029 template void
3031 const IndexSet &,
3032 const DynamicSparsityPattern &,
3033 const MPI_Comm,
3034 const bool);
3035
3036 template void
3037 SparseMatrix::vmult(MPI::Vector &, const MPI::Vector &) const;
3038
3039 template void
3041 const ::Vector<double> &) const;
3042
3043 template void
3046 const ::LinearAlgebra::distributed::Vector<double> &) const;
3047
3048# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3049# if defined(HAVE_TPETRA_INST_DOUBLE)
3050 template void
3053 const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3054# endif
3055
3056# if defined(HAVE_TPETRA_INST_FLOAT)
3057 template void
3060 const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3061# endif
3062# endif
3063
3064 template void
3067 const ::LinearAlgebra::EpetraWrappers::Vector &) const;
3068
3069 template void
3070 SparseMatrix::Tvmult(MPI::Vector &, const MPI::Vector &) const;
3071
3072 template void
3074 const ::Vector<double> &) const;
3075
3076 template void
3079 const ::LinearAlgebra::distributed::Vector<double> &) const;
3080
3081# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3082# if defined(HAVE_TPETRA_INST_DOUBLE)
3083 template void
3086 const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3087# endif
3088
3089# if defined(HAVE_TPETRA_INST_FLOAT)
3090 template void
3093 const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3094# endif
3095# endif
3096
3097 template void
3100 const ::LinearAlgebra::EpetraWrappers::Vector &) const;
3101
3102 template void
3103 SparseMatrix::vmult_add(MPI::Vector &, const MPI::Vector &) const;
3104
3105 template void
3107 const ::Vector<double> &) const;
3108
3109 template void
3112 const ::LinearAlgebra::distributed::Vector<double> &) const;
3113
3114# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3115# if defined(HAVE_TPETRA_INST_DOUBLE)
3116 template void
3119 const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3120# endif
3121
3122# if defined(HAVE_TPETRA_INST_FLOAT)
3123 template void
3126 const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3127# endif
3128# endif
3129
3130 template void
3133 const ::LinearAlgebra::EpetraWrappers::Vector &) const;
3134
3135 template void
3136 SparseMatrix::Tvmult_add(MPI::Vector &, const MPI::Vector &) const;
3137
3138 template void
3140 const ::Vector<double> &) const;
3141
3142 template void
3145 const ::LinearAlgebra::distributed::Vector<double> &) const;
3146
3147# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3148# if defined(HAVE_TPETRA_INST_DOUBLE)
3149 template void
3152 const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3153# endif
3154
3155# if defined(HAVE_TPETRA_INST_FLOAT)
3156 template void
3159 const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3160# endif
3161# endif
3162
3163 template void
3166 const ::LinearAlgebra::EpetraWrappers::Vector &) const;
3167} // namespace TrilinosWrappers
3168# endif // DOXYGEN
3169
3171
3172#endif // DEAL_II_WITH_TRILINOS
const IndexSet & row_index_set() const
size_type row_length(const size_type row) const
size_type column_number(const size_type row, const size_type index) const
size_type size() const
Definition index_set.h:1765
bool is_element(const size_type index) const
Definition index_set.h:1883
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
const Epetra_FEVector & trilinos_vector() const
const Tpetra::Vector< Number, int, types::signed_global_dof_index > & trilinos_vector() const
size_type n_rows() const
size_type n_cols() const
void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)
const Epetra_BlockMap & trilinos_partitioner() const
size_type size() const override
const Epetra_MultiVector & trilinos_vector() const
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
std::pair< size_type, size_type > local_range() const
std::shared_ptr< std::vector< TrilinosScalar > > value_cache
std::shared_ptr< std::vector< size_type > > colnum_cache
void set(const size_type i, const size_type j, const TrilinosScalar value)
std::unique_ptr< Epetra_Map > column_space_map
SparseMatrix & operator*=(const TrilinosScalar factor)
void mmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
std::unique_ptr< Epetra_Export > nonlocal_matrix_exporter
void compress(VectorOperation::values operation)
std::unique_ptr< Epetra_FECrsMatrix > matrix
const Epetra_CrsMatrix & trilinos_matrix() const
TrilinosScalar matrix_norm_square(const MPI::Vector &v) const
IndexSet locally_owned_range_indices() const
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
void Tmmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
std::enable_if_t< std::is_same_v< typename VectorType::value_type, TrilinosScalar > > Tvmult(VectorType &dst, const VectorType &src) const
void clear_row(const size_type row, const TrilinosScalar new_diag_value=0)
void clear_rows(const ArrayView< const size_type > &rows, const TrilinosScalar new_diag_value=0)
void reinit(const SparsityPatternType &sparsity_pattern)
void vmult_add(VectorType &dst, const VectorType &src) const
SparseMatrix & operator=(const SparseMatrix &)=delete
SparseMatrix & operator/=(const TrilinosScalar factor)
void Tvmult_add(VectorType &dst, const VectorType &src) const
TrilinosScalar el(const size_type i, const size_type j) const
IndexSet locally_owned_domain_indices() const
bool in_local_range(const size_type index) const
void copy_from(const SparseMatrix &source)
unsigned int row_length(const size_type row) const
std::uint64_t n_nonzero_elements() const
std::enable_if_t< std::is_same_v< typename VectorType::value_type, TrilinosScalar > > vmult(VectorType &dst, const VectorType &src) const
const Epetra_CrsGraph & trilinos_sparsity_pattern() const
TrilinosScalar diag_element(const size_type i) const
void add(const size_type i, const size_type j, const TrilinosScalar value)
unsigned int local_size() const
TrilinosScalar operator()(const size_type i, const size_type j) const
TrilinosScalar matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const
std::pair< size_type, size_type > local_range() const
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
const Epetra_FECrsGraph & trilinos_sparsity_pattern() const
::types::global_dof_index size_type
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
virtual int Apply(const VectorType &X, VectorType &Y) const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define AssertIsFinite(number)
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMatrixNotCompressed()
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1193
spacedim const Point< spacedim > & p
Definition grid_tools.h:990
@ matrix
Contents is actually a matrix.
TrilinosPayload operator+(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
TrilinosPayload operator*(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
void check_vector_map_equality(const Epetra_CrsMatrix &, const VectorType &, const VectorType &)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
void perform_mmult(const SparseMatrix &inputleft, const SparseMatrix &inputright, SparseMatrix &result, const MPI::Vector &V, const bool transpose_left)
TrilinosWrappers::types::int_type min_my_gid(const Epetra_BlockMap &map)
TrilinosWrappers::types::int64_type n_global_elements(const Epetra_BlockMap &map)
TrilinosWrappers::types::int_type global_column_index(const Epetra_CrsMatrix &matrix, const ::types::global_dof_index i)
TrilinosWrappers::types::int_type max_my_gid(const Epetra_BlockMap &map)
TrilinosWrappers::types::int_type global_row_index(const Epetra_CrsMatrix &matrix, const ::types::global_dof_index i)
TrilinosWrappers::types::int_type n_global_cols(const Epetra_CrsGraph &graph)
const Epetra_Comm & comm_self()
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1032
Definition types.h:32
double TrilinosScalar
Definition types.h:178