Reference documentation for deal.II version GIT relicensing-136-gb80d0be4af 2024-03-18 08:20: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 int ierr;
1603 if (last_action == Insert)
1604 {
1605 // TODO: this could lead to a dead lock when only one processor
1606 // calls GlobalAssemble.
1607 ierr =
1608 matrix->GlobalAssemble(*column_space_map, matrix->RowMap(), false);
1609
1610 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1611 }
1612
1613 last_action = Add;
1614
1615 const TrilinosWrappers::types::int_type *col_index_ptr;
1616 const TrilinosScalar *col_value_ptr;
1617 const TrilinosWrappers::types::int_type trilinos_row = row;
1619
1620 boost::container::small_vector<TrilinosScalar, 100> local_value_array(
1621 n_cols);
1622 boost::container::small_vector<TrilinosWrappers::types::int_type, 100>
1623 local_index_array(n_cols);
1624
1625 // If we don't elide zeros, the pointers are already available... need to
1626 // cast to non-const pointers as that is the format taken by Trilinos (but
1627 // we will not modify const data)
1628 if (elide_zero_values == false)
1629 {
1630 col_index_ptr =
1631 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1632 col_indices);
1633 col_value_ptr = values;
1634 n_columns = n_cols;
1635# ifdef DEBUG
1636 for (size_type j = 0; j < n_cols; ++j)
1637 AssertIsFinite(values[j]);
1638# endif
1639 }
1640 else
1641 {
1642 // Otherwise, extract nonzero values in each row and the corresponding
1643 // index.
1644 col_index_ptr = local_index_array.data();
1645 col_value_ptr = local_value_array.data();
1646
1647 n_columns = 0;
1648 for (size_type j = 0; j < n_cols; ++j)
1649 {
1650 const double value = values[j];
1651
1652 AssertIsFinite(value);
1653 if (value != 0)
1654 {
1655 local_index_array[n_columns] = col_indices[j];
1656 local_value_array[n_columns] = value;
1657 ++n_columns;
1658 }
1659 }
1660
1661 AssertIndexRange(n_columns, n_cols + 1);
1662 }
1663 // Exit early if there is nothing to do
1664 if (n_columns == 0)
1665 {
1666 return;
1667 }
1668
1669 // If the calling processor owns the row to which we want to add values, we
1670 // can directly call the Epetra_CrsMatrix input function, which is much
1671 // faster than the Epetra_FECrsMatrix function.
1672 if (matrix->RowMap().MyGID(
1673 static_cast<TrilinosWrappers::types::int_type>(row)) == true)
1674 {
1675 ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row,
1676 n_columns,
1677 col_value_ptr,
1678 col_index_ptr);
1679 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1680 }
1681 else if (nonlocal_matrix.get() != nullptr)
1682 {
1683 compressed = false;
1684 // this is the case when we have explicitly set the off-processor rows
1685 // and want to create a separate matrix object for them (to retain
1686 // thread-safety)
1687 Assert(nonlocal_matrix->RowMap().LID(
1688 static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1689 ExcMessage("Attempted to write into off-processor matrix row "
1690 "that has not be specified as being writable upon "
1691 "initialization"));
1692 ierr = nonlocal_matrix->SumIntoGlobalValues(row,
1693 n_columns,
1694 col_value_ptr,
1695 col_index_ptr);
1696 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1697 }
1698 else
1699 {
1700 // When we're at off-processor data, we have to stick with the
1701 // standard SumIntoGlobalValues function. Nevertheless, the way we
1702 // call it is the fastest one (any other will lead to repeated
1703 // allocation and deallocation of memory in order to call the function
1704 // we already use, which is very inefficient if writing one element at
1705 // a time).
1706 compressed = false;
1707
1708 ierr = matrix->SumIntoGlobalValues(1,
1709 &trilinos_row,
1710 n_columns,
1711 col_index_ptr,
1712 &col_value_ptr,
1713 Epetra_FECrsMatrix::ROW_MAJOR);
1714 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1715 }
1716
1717# ifdef DEBUG
1718 if (ierr > 0)
1719 {
1720 std::cout << "------------------------------------------" << std::endl;
1721 std::cout << "Got error " << ierr << " in row " << row << " of proc "
1722 << matrix->RowMap().Comm().MyPID()
1723 << " when trying to add the columns:" << std::endl;
1724 for (TrilinosWrappers::types::int_type i = 0; i < n_columns; ++i)
1725 std::cout << col_index_ptr[i] << " ";
1726 std::cout << std::endl << std::endl;
1727 std::cout << "Matrix row "
1728 << (matrix->RowMap().MyGID(
1729 static_cast<TrilinosWrappers::types::int_type>(row)) ==
1730 false ?
1731 "(nonlocal part)" :
1732 "")
1733 << " has the following indices:" << std::endl;
1734 std::vector<TrilinosWrappers::types::int_type> indices;
1735 const Epetra_CrsGraph *graph =
1736 (nonlocal_matrix.get() != nullptr &&
1737 matrix->RowMap().MyGID(
1738 static_cast<TrilinosWrappers::types::int_type>(row)) == false) ?
1739 &nonlocal_matrix->Graph() :
1740 &matrix->Graph();
1741
1742 indices.resize(graph->NumGlobalIndices(row));
1743 int n_indices = 0;
1744 graph->ExtractGlobalRowCopy(row,
1745 indices.size(),
1746 n_indices,
1747 indices.data());
1748 AssertDimension(n_indices, indices.size());
1749
1750 for (TrilinosWrappers::types::int_type i = 0; i < n_indices; ++i)
1751 std::cout << indices[i] << " ";
1752 std::cout << std::endl << std::endl;
1753 Assert(ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
1754 }
1755# endif
1756 Assert(ierr >= 0, ExcTrilinosError(ierr));
1757 }
1758
1759
1760
1761 SparseMatrix &
1763 {
1764 (void)d;
1766 compress(VectorOperation::unknown); // TODO: why do we do this? Should we
1767 // not check for is_compressed?
1768
1769 // As checked above, we are only allowed to use d==0.0, so pass
1770 // a constant zero (instead of a run-time value 'd' that *happens* to
1771 // have a zero value) to the underlying class in hopes that the compiler
1772 // can optimize this somehow.
1773 const int ierr = matrix->PutScalar(/*d=*/0.0);
1774 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1775
1776 if (nonlocal_matrix.get() != nullptr)
1777 {
1778 const int ierr = nonlocal_matrix->PutScalar(/*d=*/0.0);
1779 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1780 }
1781
1782 return *this;
1783 }
1784
1785
1786
1787 void
1789 {
1790 AssertDimension(rhs.m(), m());
1791 AssertDimension(rhs.n(), n());
1794 Assert(matrix->RowMap().SameAs(rhs.matrix->RowMap()),
1795 ExcMessage("Can only add matrices with same distribution of rows"));
1796 Assert(matrix->Filled() && rhs.matrix->Filled(),
1797 ExcMessage("Addition of matrices only allowed if matrices are "
1798 "filled, i.e., compress() has been called"));
1799
1800 const bool same_col_map = matrix->ColMap().SameAs(rhs.matrix->ColMap());
1801
1802 for (const auto row : locally_owned_range_indices())
1803 {
1804 const int row_local = matrix->RowMap().LID(
1805 static_cast<TrilinosWrappers::types::int_type>(row));
1806 Assert((row_local >= 0), ExcAccessToNonlocalRow(row));
1807
1808 // First get a view to the matrix columns of both matrices. Note that
1809 // the data is in local index spaces so we need to be careful not only
1810 // to compare column indices in case they are derived from the same
1811 // map.
1812 int n_entries, rhs_n_entries;
1813 TrilinosScalar *value_ptr, *rhs_value_ptr;
1814 int *index_ptr, *rhs_index_ptr;
1815 int ierr = rhs.matrix->ExtractMyRowView(row_local,
1816 rhs_n_entries,
1817 rhs_value_ptr,
1818 rhs_index_ptr);
1819 (void)ierr;
1820 Assert(ierr == 0, ExcTrilinosError(ierr));
1821
1822 ierr =
1823 matrix->ExtractMyRowView(row_local, n_entries, value_ptr, index_ptr);
1824 Assert(ierr == 0, ExcTrilinosError(ierr));
1825 bool expensive_checks = (n_entries != rhs_n_entries || !same_col_map);
1826 if (!expensive_checks)
1827 {
1828 // check if the column indices are the same. If yes, can simply
1829 // copy over the data.
1830 expensive_checks = std::memcmp(static_cast<void *>(index_ptr),
1831 static_cast<void *>(rhs_index_ptr),
1832 sizeof(int) * n_entries) != 0;
1833 if (!expensive_checks)
1834 for (int i = 0; i < n_entries; ++i)
1835 value_ptr[i] += rhs_value_ptr[i] * factor;
1836 }
1837 // Now to the expensive case where we need to check all column indices
1838 // against each other (transformed into global index space) and where
1839 // we need to make sure that all entries we are about to add into the
1840 // lhs matrix actually exist
1841 if (expensive_checks)
1842 {
1843 for (int i = 0; i < rhs_n_entries; ++i)
1844 {
1845 if (rhs_value_ptr[i] == 0.)
1846 continue;
1847 const TrilinosWrappers::types::int_type rhs_global_col =
1848 global_column_index(*rhs.matrix, rhs_index_ptr[i]);
1849 int local_col = matrix->ColMap().LID(rhs_global_col);
1850 int *local_index = Utilities::lower_bound(index_ptr,
1851 index_ptr + n_entries,
1852 local_col);
1853 Assert(local_index != index_ptr + n_entries &&
1854 *local_index == local_col,
1855 ExcMessage(
1856 "Adding the entries from the other matrix "
1857 "failed, because the sparsity pattern "
1858 "of that matrix includes more elements than the "
1859 "calling matrix, which is not allowed."));
1860 value_ptr[local_index - index_ptr] += factor * rhs_value_ptr[i];
1861 }
1862 }
1863 }
1864 }
1865
1866
1867
1868 void
1870 {
1871 // This only flips a flag that tells
1872 // Trilinos that any vmult operation
1873 // should be done with the
1874 // transpose. However, the matrix
1875 // structure is not reset.
1876 int ierr;
1877
1878 if (!matrix->UseTranspose())
1879 {
1880 ierr = matrix->SetUseTranspose(true);
1881 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1882 }
1883 else
1884 {
1885 ierr = matrix->SetUseTranspose(false);
1886 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1887 }
1888 }
1889
1890
1891
1892 SparseMatrix &
1894 {
1895 const int ierr = matrix->Scale(a);
1896 Assert(ierr == 0, ExcTrilinosError(ierr));
1897 (void)ierr; // removes -Wunused-variable in optimized mode
1898
1899 return *this;
1900 }
1901
1902
1903
1904 SparseMatrix &
1906 {
1907 Assert(a != 0, ExcDivideByZero());
1908
1909 const TrilinosScalar factor = 1. / a;
1910
1911 const int ierr = matrix->Scale(factor);
1912 Assert(ierr == 0, ExcTrilinosError(ierr));
1913 (void)ierr; // removes -Wunused-variable in optimized mode
1914
1915 return *this;
1916 }
1917
1918
1919
1922 {
1923 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1924 return matrix->NormOne();
1925 }
1926
1927
1928
1931 {
1932 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1933 return matrix->NormInf();
1934 }
1935
1936
1937
1940 {
1941 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1942 return matrix->NormFrobenius();
1943 }
1944
1945
1946
1947 namespace internal
1948 {
1949 namespace SparseMatrixImplementation
1950 {
1951 template <typename VectorType>
1952 inline void
1953 check_vector_map_equality(const Epetra_CrsMatrix &,
1954 const VectorType &,
1955 const VectorType &)
1956 {}
1957
1958 inline void
1959 check_vector_map_equality(const Epetra_CrsMatrix &m,
1962 {
1963 Assert(in.trilinos_partitioner().SameAs(m.DomainMap()) == true,
1964 ExcMessage("The column partitioning of a matrix does not match "
1965 "the partitioning of a vector you are trying to "
1966 "multiply it with. Are you multiplying the "
1967 "matrix with a vector that has ghost elements?"));
1968 Assert(out.trilinos_partitioner().SameAs(m.RangeMap()) == true,
1969 ExcMessage("The row partitioning of a matrix does not match "
1970 "the partitioning of a vector you are trying to "
1971 "put the result of a matrix-vector product in. "
1972 "Are you trying to put the product of the "
1973 "matrix with a vector into a vector that has "
1974 "ghost elements?"));
1975 (void)m;
1976 (void)in;
1977 (void)out;
1978 }
1979 } // namespace SparseMatrixImplementation
1980 } // namespace internal
1981
1982
1983 template <typename VectorType>
1984 std::enable_if_t<
1985 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1986 SparseMatrix::vmult(VectorType &dst, const VectorType &src) const
1987 {
1988 Assert(&src != &dst, ExcSourceEqualsDestination());
1989 Assert(matrix->Filled(), ExcMatrixNotCompressed());
1990 (void)src;
1991 (void)dst;
1992
1994 src,
1995 dst);
1996 const size_type dst_local_size = internal::end(dst) - internal::begin(dst);
1997 AssertDimension(dst_local_size, matrix->RangeMap().NumMyPoints());
1998 const size_type src_local_size = internal::end(src) - internal::begin(src);
1999 AssertDimension(src_local_size, matrix->DomainMap().NumMyPoints());
2000
2001 Epetra_MultiVector tril_dst(
2002 View, matrix->RangeMap(), internal::begin(dst), dst_local_size, 1);
2003 Epetra_MultiVector tril_src(View,
2004 matrix->DomainMap(),
2005 const_cast<TrilinosScalar *>(
2006 internal::begin(src)),
2007 src_local_size,
2008 1);
2009
2010 const int ierr = matrix->Multiply(false, tril_src, tril_dst);
2011 Assert(ierr == 0, ExcTrilinosError(ierr));
2012 (void)ierr; // removes -Wunused-variable in optimized mode
2013 }
2014
2015
2016
2017 template <typename VectorType>
2018 std::enable_if_t<
2019 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2020 SparseMatrix::vmult(VectorType & /*dst*/, const VectorType & /*src*/) const
2021 {
2023 }
2024
2025
2026
2027 template <typename VectorType>
2028 std::enable_if_t<
2029 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2030 SparseMatrix::Tvmult(VectorType &dst, const VectorType &src) const
2031 {
2032 Assert(&src != &dst, ExcSourceEqualsDestination());
2033 Assert(matrix->Filled(), ExcMatrixNotCompressed());
2034
2036 dst,
2037 src);
2038 const size_type dst_local_size = internal::end(dst) - internal::begin(dst);
2039 AssertDimension(dst_local_size, matrix->DomainMap().NumMyPoints());
2040 const size_type src_local_size = internal::end(src) - internal::begin(src);
2041 AssertDimension(src_local_size, matrix->RangeMap().NumMyPoints());
2042
2043 Epetra_MultiVector tril_dst(
2044 View, matrix->DomainMap(), internal::begin(dst), dst_local_size, 1);
2045 Epetra_MultiVector tril_src(View,
2046 matrix->RangeMap(),
2047 const_cast<double *>(internal::begin(src)),
2048 src_local_size,
2049 1);
2050
2051 const int ierr = matrix->Multiply(true, tril_src, tril_dst);
2052 Assert(ierr == 0, ExcTrilinosError(ierr));
2053 (void)ierr; // removes -Wunused-variable in optimized mode
2054 }
2055
2056
2057
2058 template <typename VectorType>
2059 std::enable_if_t<
2060 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2061 SparseMatrix::Tvmult(VectorType & /*dst*/, const VectorType & /*src*/) const
2062 {
2064 }
2065
2066
2067
2068 template <typename VectorType>
2069 void
2070 SparseMatrix::vmult_add(VectorType &dst, const VectorType &src) const
2071 {
2072 Assert(&src != &dst, ExcSourceEqualsDestination());
2073
2074 // Reinit a temporary vector with fast argument set, which does not
2075 // overwrite the content (to save time).
2076 VectorType tmp_vector;
2077 tmp_vector.reinit(dst, true);
2078 vmult(tmp_vector, src);
2079 dst += tmp_vector;
2080 }
2081
2082
2083
2084 template <typename VectorType>
2085 void
2086 SparseMatrix::Tvmult_add(VectorType &dst, const VectorType &src) const
2087 {
2088 Assert(&src != &dst, ExcSourceEqualsDestination());
2089
2090 // Reinit a temporary vector with fast argument set, which does not
2091 // overwrite the content (to save time).
2092 VectorType tmp_vector;
2093 tmp_vector.reinit(dst, true);
2094 Tvmult(tmp_vector, src);
2095 dst += tmp_vector;
2096 }
2097
2098
2099
2102 {
2103 AssertDimension(m(), v.size());
2104 Assert(matrix->RowMap().SameAs(matrix->DomainMap()), ExcNotQuadratic());
2105
2106 MPI::Vector temp_vector;
2107 temp_vector.reinit(v, true);
2108
2109 vmult(temp_vector, v);
2110 return temp_vector * v;
2111 }
2112
2113
2114
2117 const MPI::Vector &v) const
2118 {
2119 AssertDimension(m(), u.size());
2120 AssertDimension(m(), v.size());
2121 Assert(matrix->RowMap().SameAs(matrix->DomainMap()), ExcNotQuadratic());
2122
2123 MPI::Vector temp_vector;
2124 temp_vector.reinit(v, true);
2125
2126 vmult(temp_vector, v);
2127 return u * temp_vector;
2128 }
2129
2130
2131
2132 namespace internals
2133 {
2135
2136 void
2137 perform_mmult(const SparseMatrix &inputleft,
2138 const SparseMatrix &inputright,
2139 SparseMatrix &result,
2140 const MPI::Vector &V,
2141 const bool transpose_left)
2142 {
2143 const bool use_vector = (V.size() == inputright.m() ? true : false);
2144 if (transpose_left == false)
2145 {
2146 Assert(inputleft.n() == inputright.m(),
2147 ExcDimensionMismatch(inputleft.n(), inputright.m()));
2148 Assert(inputleft.trilinos_matrix().DomainMap().SameAs(
2149 inputright.trilinos_matrix().RangeMap()),
2150 ExcMessage("Parallel partitioning of A and B does not fit."));
2151 }
2152 else
2153 {
2154 Assert(inputleft.m() == inputright.m(),
2155 ExcDimensionMismatch(inputleft.m(), inputright.m()));
2156 Assert(inputleft.trilinos_matrix().RangeMap().SameAs(
2157 inputright.trilinos_matrix().RangeMap()),
2158 ExcMessage("Parallel partitioning of A and B does not fit."));
2159 }
2160
2161 result.clear();
2162
2163 // create a suitable operator B: in case
2164 // we do not use a vector, all we need to
2165 // do is to set the pointer. Otherwise,
2166 // we insert the data from B, but
2167 // multiply each row with the respective
2168 // vector element.
2169 Teuchos::RCP<Epetra_CrsMatrix> mod_B;
2170 if (use_vector == false)
2171 {
2172 mod_B = Teuchos::rcp(const_cast<Epetra_CrsMatrix *>(
2173 &inputright.trilinos_matrix()),
2174 false);
2175 }
2176 else
2177 {
2178 mod_B = Teuchos::rcp(
2179 new Epetra_CrsMatrix(Copy, inputright.trilinos_sparsity_pattern()),
2180 true);
2181 mod_B->FillComplete(inputright.trilinos_matrix().DomainMap(),
2182 inputright.trilinos_matrix().RangeMap());
2183 Assert(inputright.local_range() == V.local_range(),
2184 ExcMessage("Parallel distribution of matrix B and vector V "
2185 "does not match."));
2186
2187 const int local_N = inputright.local_size();
2188 for (int i = 0; i < local_N; ++i)
2189 {
2190 int N_entries = -1;
2191 double *new_data, *B_data;
2192 mod_B->ExtractMyRowView(i, N_entries, new_data);
2193 inputright.trilinos_matrix().ExtractMyRowView(i,
2194 N_entries,
2195 B_data);
2196 double value = V.trilinos_vector()[0][i];
2197 for (TrilinosWrappers::types::int_type j = 0; j < N_entries; ++j)
2198 new_data[j] = value * B_data[j];
2199 }
2200 }
2201
2202
2203 SparseMatrix tmp_result(transpose_left ?
2204 inputleft.locally_owned_domain_indices() :
2205 inputleft.locally_owned_range_indices(),
2206 inputright.locally_owned_domain_indices(),
2207 inputleft.get_mpi_communicator());
2208
2209# ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT
2210 EpetraExt::MatrixMatrix::Multiply(inputleft.trilinos_matrix(),
2211 transpose_left,
2212 *mod_B,
2213 false,
2214 const_cast<Epetra_CrsMatrix &>(
2215 tmp_result.trilinos_matrix()));
2216# else
2217 Assert(false,
2218 ExcMessage("This function requires that the Trilinos "
2219 "installation found while running the deal.II "
2220 "CMake scripts contains the optional Trilinos "
2221 "package 'EpetraExt'. However, this optional "
2222 "part of Trilinos was not found."));
2223# endif
2224 result.reinit(tmp_result.trilinos_matrix());
2225 }
2226 } // namespace internals
2227
2228
2229 void
2231 const SparseMatrix &B,
2232 const MPI::Vector &V) const
2233 {
2234 internals::perform_mmult(*this, B, C, V, false);
2235 }
2236
2237
2238
2239 void
2241 const SparseMatrix &B,
2242 const MPI::Vector &V) const
2243 {
2244 internals::perform_mmult(*this, B, C, V, true);
2245 }
2246
2247
2248
2249 void
2254
2255
2256
2257 // As of now, no particularly neat
2258 // output is generated in case of
2259 // multiple processors.
2260 void
2261 SparseMatrix::print(std::ostream &out,
2262 const bool print_detailed_trilinos_information) const
2263 {
2264 if (print_detailed_trilinos_information == true)
2265 out << *matrix;
2266 else
2267 {
2268 double *values;
2269 int *indices;
2270 int num_entries;
2271
2272 for (int i = 0; i < matrix->NumMyRows(); ++i)
2273 {
2274 const int ierr =
2275 matrix->ExtractMyRowView(i, num_entries, values, indices);
2276 (void)ierr;
2277 Assert(ierr == 0, ExcTrilinosError(ierr));
2278
2279 for (TrilinosWrappers::types::int_type j = 0; j < num_entries; ++j)
2281 << ","
2283 << ") " << values[j] << std::endl;
2284 }
2285 }
2286
2287 AssertThrow(out.fail() == false, ExcIO());
2288 }
2289
2290
2291
2294 {
2295 size_type static_memory =
2296 sizeof(*this) + sizeof(*matrix) + sizeof(*matrix->Graph().DataPtr());
2297 return (
2299 matrix->NumMyNonzeros() +
2300 sizeof(int) * local_size() + static_memory);
2301 }
2302
2303
2304
2305 MPI_Comm
2307 {
2308 const Epetra_MpiComm *mpi_comm =
2309 dynamic_cast<const Epetra_MpiComm *>(&matrix->RangeMap().Comm());
2310 Assert(mpi_comm != nullptr, ExcInternalError());
2311 return mpi_comm->Comm();
2312 }
2313} // namespace TrilinosWrappers
2314
2315
2316namespace TrilinosWrappers
2317{
2318 namespace internal
2319 {
2320 namespace LinearOperatorImplementation
2321 {
2323 : use_transpose(false)
2324 , communicator(MPI_COMM_SELF)
2325 , domain_map(IndexSet().make_trilinos_map(communicator.Comm()))
2326 , range_map(IndexSet().make_trilinos_map(communicator.Comm()))
2327 {
2328 vmult = [](Range &, const Domain &) {
2329 Assert(false,
2330 ExcMessage("Uninitialized TrilinosPayload::vmult called "
2331 "(Default constructor)"));
2332 };
2333
2334 Tvmult = [](Domain &, const Range &) {
2335 Assert(false,
2336 ExcMessage("Uninitialized TrilinosPayload::Tvmult called "
2337 "(Default constructor)"));
2338 };
2339
2340 inv_vmult = [](Domain &, const Range &) {
2341 Assert(false,
2342 ExcMessage("Uninitialized TrilinosPayload::inv_vmult called "
2343 "(Default constructor)"));
2344 };
2345
2346 inv_Tvmult = [](Range &, const Domain &) {
2347 Assert(false,
2348 ExcMessage("Uninitialized TrilinosPayload::inv_Tvmult called "
2349 "(Default constructor)"));
2350 };
2351 }
2352
2353
2354
2356 const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2357 const TrilinosWrappers::SparseMatrix &matrix)
2358 : TrilinosPayload(const_cast<Epetra_CrsMatrix &>(
2359 matrix.trilinos_matrix()),
2360 /*op_supports_inverse_operations = */ false,
2361 matrix_exemplar.trilinos_matrix().UseTranspose(),
2362 matrix_exemplar.get_mpi_communicator(),
2363 matrix_exemplar.locally_owned_domain_indices(),
2364 matrix_exemplar.locally_owned_range_indices())
2365 {}
2366
2367
2368
2370 const TrilinosPayload &payload_exemplar,
2371 const TrilinosWrappers::SparseMatrix &matrix)
2372
2373 : TrilinosPayload(const_cast<Epetra_CrsMatrix &>(
2374 matrix.trilinos_matrix()),
2375 /*op_supports_inverse_operations = */ false,
2376 payload_exemplar.UseTranspose(),
2377 payload_exemplar.get_mpi_communicator(),
2378 payload_exemplar.locally_owned_domain_indices(),
2379 payload_exemplar.locally_owned_range_indices())
2380 {}
2381
2382
2383
2385 const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2386 const TrilinosWrappers::PreconditionBase &preconditioner)
2387 : TrilinosPayload(preconditioner.trilinos_operator(),
2388 /*op_supports_inverse_operations = */ true,
2389 matrix_exemplar.trilinos_matrix().UseTranspose(),
2390 matrix_exemplar.get_mpi_communicator(),
2391 matrix_exemplar.locally_owned_domain_indices(),
2392 matrix_exemplar.locally_owned_range_indices())
2393 {}
2394
2395
2396
2398 const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2399 const TrilinosWrappers::PreconditionBase &preconditioner)
2401 preconditioner.trilinos_operator(),
2402 /*op_supports_inverse_operations = */ true,
2403 preconditioner_exemplar.trilinos_operator().UseTranspose(),
2404 preconditioner_exemplar.get_mpi_communicator(),
2405 preconditioner_exemplar.locally_owned_domain_indices(),
2406 preconditioner_exemplar.locally_owned_range_indices())
2407 {}
2408
2409
2410
2412 const TrilinosPayload &payload_exemplar,
2413 const TrilinosWrappers::PreconditionBase &preconditioner)
2414 : TrilinosPayload(preconditioner.trilinos_operator(),
2415 /*op_supports_inverse_operations = */ true,
2416 payload_exemplar.UseTranspose(),
2417 payload_exemplar.get_mpi_communicator(),
2418 payload_exemplar.locally_owned_domain_indices(),
2419 payload_exemplar.locally_owned_range_indices())
2420 {}
2421
2422
2423
2425 : vmult(payload.vmult)
2426 , Tvmult(payload.Tvmult)
2427 , inv_vmult(payload.inv_vmult)
2428 , inv_Tvmult(payload.inv_Tvmult)
2429 , use_transpose(payload.use_transpose)
2430 , communicator(payload.communicator)
2431 , domain_map(payload.domain_map)
2432 , range_map(payload.range_map)
2433 {}
2434
2435
2436
2437 // Composite copy constructor
2438 // This is required for PackagedOperations
2440 const TrilinosPayload &second_op)
2441 : use_transpose(false)
2442 , // The combination of operators provides the exact
2443 // definition of the operation
2444 communicator(first_op.communicator)
2445 , domain_map(second_op.domain_map)
2446 , range_map(first_op.range_map)
2447 {}
2448
2449
2450
2453 {
2454 TrilinosPayload return_op(*this);
2455
2456 return_op.vmult = [](Range &tril_dst, const Range &tril_src) {
2457 tril_dst = tril_src;
2458 };
2459
2460 return_op.Tvmult = [](Range &tril_dst, const Range &tril_src) {
2461 tril_dst = tril_src;
2462 };
2463
2464 return_op.inv_vmult = [](Range &tril_dst, const Range &tril_src) {
2465 tril_dst = tril_src;
2466 };
2467
2468 return_op.inv_Tvmult = [](Range &tril_dst, const Range &tril_src) {
2469 tril_dst = tril_src;
2470 };
2471
2472 return return_op;
2473 }
2474
2475
2476
2479 {
2480 TrilinosPayload return_op(*this);
2481
2482 return_op.vmult = [](Range &tril_dst, const Domain &) {
2483 const int ierr = tril_dst.PutScalar(0.0);
2484
2485 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2486 };
2487
2488 return_op.Tvmult = [](Domain &tril_dst, const Range &) {
2489 const int ierr = tril_dst.PutScalar(0.0);
2490
2491 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2492 };
2493
2494 return_op.inv_vmult = [](Domain &tril_dst, const Range &) {
2495 AssertThrow(false,
2496 ExcMessage("Cannot compute inverse of null operator"));
2497
2498 const int ierr = tril_dst.PutScalar(0.0);
2499
2500 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2501 };
2502
2503 return_op.inv_Tvmult = [](Range &tril_dst, const Domain &) {
2504 AssertThrow(false,
2505 ExcMessage("Cannot compute inverse of null operator"));
2506
2507 const int ierr = tril_dst.PutScalar(0.0);
2508
2509 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2510 };
2511
2512 return return_op;
2513 }
2514
2515
2516
2519 {
2520 TrilinosPayload return_op(*this);
2521 return_op.transpose();
2522 return return_op;
2523 }
2524
2525
2526
2527 IndexSet
2532
2533
2534
2535 IndexSet
2540
2541
2542
2543 MPI_Comm
2545 {
2546 return communicator.Comm();
2547 }
2548
2549
2550
2551 void
2556
2557
2558
2559 bool
2561 {
2562 return use_transpose;
2563 }
2564
2565
2566
2567 int
2569 {
2571 {
2573 std::swap(domain_map, range_map);
2574 std::swap(vmult, Tvmult);
2575 std::swap(inv_vmult, inv_Tvmult);
2576 }
2577 return 0;
2578 }
2579
2580
2581
2582 int
2584 {
2585 // The transposedness of the operations is taken care of
2586 // when we hit the transpose flag.
2587 vmult(Y, X);
2588 return 0;
2589 }
2590
2591
2592
2593 int
2595 {
2596 // The transposedness of the operations is taken care of
2597 // when we hit the transpose flag.
2598 inv_vmult(X, Y);
2599 return 0;
2600 }
2601
2602
2603
2604 const char *
2606 {
2607 return "TrilinosPayload";
2608 }
2609
2610
2611
2612 const Epetra_Comm &
2614 {
2615 return communicator;
2616 }
2617
2618
2619
2620 const Epetra_Map &
2622 {
2623 return domain_map;
2624 }
2625
2626
2627
2628 const Epetra_Map &
2630 {
2631 return range_map;
2632 }
2633
2634
2635
2636 bool
2638 {
2639 return false;
2640 }
2641
2642
2643
2644 double
2646 {
2648 return 0.0;
2649 }
2650
2651
2652
2655 const TrilinosPayload &second_op)
2656 {
2657 using Domain = typename TrilinosPayload::Domain;
2658 using Range = typename TrilinosPayload::Range;
2659 using Intermediate = typename TrilinosPayload::VectorType;
2660 using GVMVectorType = TrilinosWrappers::MPI::Vector;
2661
2663 second_op.locally_owned_domain_indices(),
2664 ExcMessage(
2665 "Operators are set to work on incompatible IndexSets."));
2667 second_op.locally_owned_range_indices(),
2668 ExcMessage(
2669 "Operators are set to work on incompatible IndexSets."));
2670
2671 TrilinosPayload return_op(first_op, second_op);
2672
2673 // Capture by copy so the payloads are always valid
2674 return_op.vmult = [first_op, second_op](Range &tril_dst,
2675 const Domain &tril_src) {
2676 // Duplicated from LinearOperator::operator*
2677 // TODO: Template the constructor on GrowingVectorMemory vector type?
2679 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2680
2681 // Initialize intermediate vector
2682 const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2683 i->reinit(IndexSet(first_op_init_map),
2684 first_op.get_mpi_communicator(),
2685 /*bool omit_zeroing_entries =*/true);
2686
2687 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2688 const size_type i_local_size = i->end() - i->begin();
2689 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2690 const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2691 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2692 (void)second_op_init_map;
2693 Intermediate tril_int(View,
2694 first_op_init_map,
2695 const_cast<TrilinosScalar *>(i->begin()),
2696 i_local_size,
2697 1);
2698
2699 // These operators may themselves be transposed or not, so we let them
2700 // decide what the intended outcome is
2701 second_op.Apply(tril_src, tril_int);
2702 first_op.Apply(tril_src, tril_dst);
2703 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2704 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2705 };
2706
2707 return_op.Tvmult = [first_op, second_op](Domain &tril_dst,
2708 const Range &tril_src) {
2709 // Duplicated from LinearOperator::operator*
2710 // TODO: Template the constructor on GrowingVectorMemory vector type?
2712 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2713
2714 // These operators may themselves be transposed or not, so we let them
2715 // decide what the intended outcome is
2716 // We must first transpose the operators to get the right IndexSets
2717 // for the input, intermediate and result vectors
2718 const_cast<TrilinosPayload &>(first_op).transpose();
2719 const_cast<TrilinosPayload &>(second_op).transpose();
2720
2721 // Initialize intermediate vector
2722 const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2723 i->reinit(IndexSet(first_op_init_map),
2724 first_op.get_mpi_communicator(),
2725 /*bool omit_zeroing_entries =*/true);
2726
2727 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2728 const size_type i_local_size = i->end() - i->begin();
2729 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2730 const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2731 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2732 (void)second_op_init_map;
2733 Intermediate tril_int(View,
2734 first_op_init_map,
2735 const_cast<TrilinosScalar *>(i->begin()),
2736 i_local_size,
2737 1);
2738
2739 // These operators may themselves be transposed or not, so we let them
2740 // decide what the intended outcome is
2741 second_op.Apply(tril_src, tril_int);
2742 first_op.Apply(tril_src, tril_dst);
2743 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2744 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2745
2746 // Reset transpose flag
2747 const_cast<TrilinosPayload &>(first_op).transpose();
2748 const_cast<TrilinosPayload &>(second_op).transpose();
2749 };
2750
2751 return_op.inv_vmult = [first_op, second_op](Domain &tril_dst,
2752 const Range &tril_src) {
2753 // Duplicated from LinearOperator::operator*
2754 // TODO: Template the constructor on GrowingVectorMemory vector type?
2756 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2757
2758 // Initialize intermediate vector
2759 const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2760 i->reinit(IndexSet(first_op_init_map),
2761 first_op.get_mpi_communicator(),
2762 /*bool omit_zeroing_entries =*/true);
2763
2764 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2765 const size_type i_local_size = i->end() - i->begin();
2766 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2767 const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2768 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2769 (void)second_op_init_map;
2770 Intermediate tril_int(View,
2771 first_op_init_map,
2772 const_cast<TrilinosScalar *>(i->begin()),
2773 i_local_size,
2774 1);
2775
2776 // These operators may themselves be transposed or not, so we let them
2777 // decide what the intended outcome is
2778 second_op.ApplyInverse(tril_src, tril_int);
2779 first_op.ApplyInverse(tril_src, tril_dst);
2780 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2781 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2782 };
2783
2784 return_op.inv_Tvmult = [first_op, second_op](Range &tril_dst,
2785 const Domain &tril_src) {
2786 // Duplicated from LinearOperator::operator*
2787 // TODO: Template the constructor on GrowingVectorMemory vector type?
2789 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2790
2791 // These operators may themselves be transposed or not, so we let them
2792 // decide what the intended outcome is
2793 // We must first transpose the operators to get the right IndexSets
2794 // for the input, intermediate and result vectors
2795 const_cast<TrilinosPayload &>(first_op).transpose();
2796 const_cast<TrilinosPayload &>(second_op).transpose();
2797
2798 // Initialize intermediate vector
2799 const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2800 i->reinit(IndexSet(first_op_init_map),
2801 first_op.get_mpi_communicator(),
2802 /*bool omit_zeroing_entries =*/true);
2803
2804 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2805 const size_type i_local_size = i->end() - i->begin();
2806 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2807 const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2808 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2809 (void)second_op_init_map;
2810 Intermediate tril_int(View,
2811 first_op_init_map,
2812 const_cast<TrilinosScalar *>(i->begin()),
2813 i_local_size,
2814 1);
2815
2816 // These operators may themselves be transposed or not, so we let them
2817 // decide what the intended outcome is
2818 second_op.ApplyInverse(tril_src, tril_int);
2819 first_op.ApplyInverse(tril_src, tril_dst);
2820 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2821 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2822
2823 // Reset transpose flag
2824 const_cast<TrilinosPayload &>(first_op).transpose();
2825 const_cast<TrilinosPayload &>(second_op).transpose();
2826 };
2827
2828 return return_op;
2829 }
2830
2831
2832
2833 TrilinosPayload
2835 const TrilinosPayload &second_op)
2836 {
2837 using Domain = typename TrilinosPayload::Domain;
2838 using Range = typename TrilinosPayload::Range;
2839 using Intermediate = typename TrilinosPayload::VectorType;
2840 using GVMVectorType = TrilinosWrappers::MPI::Vector;
2841
2843 second_op.locally_owned_range_indices(),
2844 ExcMessage(
2845 "Operators are set to work on incompatible IndexSets."));
2846
2847 TrilinosPayload return_op(first_op, second_op);
2848
2849 // Capture by copy so the payloads are always valid
2850 return_op.vmult = [first_op, second_op](Range &tril_dst,
2851 const Domain &tril_src) {
2852 // Duplicated from LinearOperator::operator*
2853 // TODO: Template the constructor on GrowingVectorMemory vector type?
2855 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2856
2857 // Initialize intermediate vector
2858 const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2859 i->reinit(IndexSet(first_op_init_map),
2860 first_op.get_mpi_communicator(),
2861 /*bool omit_zeroing_entries =*/true);
2862
2863 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2864 const size_type i_local_size = i->end() - i->begin();
2865 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2866 const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2867 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2868 (void)second_op_init_map;
2869 Intermediate tril_int(View,
2870 first_op_init_map,
2871 const_cast<TrilinosScalar *>(i->begin()),
2872 i_local_size,
2873 1);
2874
2875 // These operators may themselves be transposed or not, so we let them
2876 // decide what the intended outcome is
2877 second_op.Apply(tril_src, tril_int);
2878 first_op.Apply(tril_int, tril_dst);
2879 };
2880
2881 return_op.Tvmult = [first_op, second_op](Domain &tril_dst,
2882 const Range &tril_src) {
2883 // Duplicated from LinearOperator::operator*
2884 // TODO: Template the constructor on GrowingVectorMemory vector type?
2886 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2887
2888 // These operators may themselves be transposed or not, so we let them
2889 // decide what the intended outcome is
2890 // We must first transpose the operators to get the right IndexSets
2891 // for the input, intermediate and result vectors
2892 const_cast<TrilinosPayload &>(first_op).transpose();
2893 const_cast<TrilinosPayload &>(second_op).transpose();
2894
2895 // Initialize intermediate vector
2896 const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2897 i->reinit(IndexSet(first_op_init_map),
2898 first_op.get_mpi_communicator(),
2899 /*bool omit_zeroing_entries =*/true);
2900
2901 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2902 const size_type i_local_size = i->end() - i->begin();
2903 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2904 const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2905 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2906 (void)second_op_init_map;
2907 Intermediate tril_int(View,
2908 first_op_init_map,
2909 const_cast<TrilinosScalar *>(i->begin()),
2910 i_local_size,
2911 1);
2912
2913 // Apply the operators in the reverse order to vmult
2914 first_op.Apply(tril_src, tril_int);
2915 second_op.Apply(tril_int, tril_dst);
2916
2917 // Reset transpose flag
2918 const_cast<TrilinosPayload &>(first_op).transpose();
2919 const_cast<TrilinosPayload &>(second_op).transpose();
2920 };
2921
2922 return_op.inv_vmult = [first_op, second_op](Domain &tril_dst,
2923 const Range &tril_src) {
2924 // Duplicated from LinearOperator::operator*
2925 // TODO: Template the constructor on GrowingVectorMemory vector type?
2927 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2928
2929 // Initialize intermediate vector
2930 const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2931 i->reinit(IndexSet(first_op_init_map),
2932 first_op.get_mpi_communicator(),
2933 /*bool omit_zeroing_entries =*/true);
2934
2935 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2936 const size_type i_local_size = i->end() - i->begin();
2937 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2938 const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2939 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2940 (void)second_op_init_map;
2941 Intermediate tril_int(View,
2942 first_op_init_map,
2943 const_cast<TrilinosScalar *>(i->begin()),
2944 i_local_size,
2945 1);
2946
2947 // Apply the operators in the reverse order to vmult
2948 // and the same order as Tvmult
2949 first_op.ApplyInverse(tril_src, tril_int);
2950 second_op.ApplyInverse(tril_int, tril_dst);
2951 };
2952
2953 return_op.inv_Tvmult = [first_op, second_op](Range &tril_dst,
2954 const Domain &tril_src) {
2955 // Duplicated from LinearOperator::operator*
2956 // TODO: Template the constructor on GrowingVectorMemory vector type?
2958 VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2959
2960 // These operators may themselves be transposed or not, so we let them
2961 // decide what the intended outcome is
2962 // We must first transpose the operators to get the right IndexSets
2963 // for the input, intermediate and result vectors
2964 const_cast<TrilinosPayload &>(first_op).transpose();
2965 const_cast<TrilinosPayload &>(second_op).transpose();
2966
2967 // Initialize intermediate vector
2968 const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2969 i->reinit(IndexSet(first_op_init_map),
2970 first_op.get_mpi_communicator(),
2971 /*bool omit_zeroing_entries =*/true);
2972
2973 // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2974 const size_type i_local_size = i->end() - i->begin();
2975 AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2976 const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2977 AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2978 (void)second_op_init_map;
2979 Intermediate tril_int(View,
2980 first_op_init_map,
2981 const_cast<TrilinosScalar *>(i->begin()),
2982 i_local_size,
2983 1);
2984
2985 // These operators may themselves be transposed or not, so we let them
2986 // decide what the intended outcome is
2987 // Apply the operators in the reverse order to Tvmult
2988 // and the same order as vmult
2989 second_op.ApplyInverse(tril_src, tril_int);
2990 first_op.ApplyInverse(tril_int, tril_dst);
2991
2992 // Reset transpose flag
2993 const_cast<TrilinosPayload &>(first_op).transpose();
2994 const_cast<TrilinosPayload &>(second_op).transpose();
2995 };
2996
2997 return return_op;
2998 }
2999
3000 } // namespace LinearOperatorImplementation
3001 } /* namespace internal */
3002} /* namespace TrilinosWrappers */
3003
3004
3005
3006// explicit instantiations
3007# include "trilinos_sparse_matrix.inst"
3008
3009# ifndef DOXYGEN
3010// TODO: put these instantiations into generic file
3011namespace TrilinosWrappers
3012{
3013 template void
3014 SparseMatrix::reinit(const ::SparsityPattern &);
3015
3016 template void
3018
3019 template void
3021 const IndexSet &,
3022 const ::SparsityPattern &,
3023 const MPI_Comm,
3024 const bool);
3025
3026 template void
3028 const IndexSet &,
3029 const DynamicSparsityPattern &,
3030 const MPI_Comm,
3031 const bool);
3032
3033 template void
3034 SparseMatrix::vmult(MPI::Vector &, const MPI::Vector &) const;
3035
3036 template void
3038 const ::Vector<double> &) const;
3039
3040 template void
3043 const ::LinearAlgebra::distributed::Vector<double> &) const;
3044
3045# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3046# if defined(HAVE_TPETRA_INST_DOUBLE)
3047 template void
3050 const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3051# endif
3052
3053# if defined(HAVE_TPETRA_INST_FLOAT)
3054 template void
3057 const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3058# endif
3059# endif
3060
3061 template void
3064 const ::LinearAlgebra::EpetraWrappers::Vector &) const;
3065
3066 template void
3067 SparseMatrix::Tvmult(MPI::Vector &, const MPI::Vector &) const;
3068
3069 template void
3071 const ::Vector<double> &) const;
3072
3073 template void
3076 const ::LinearAlgebra::distributed::Vector<double> &) const;
3077
3078# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3079# if defined(HAVE_TPETRA_INST_DOUBLE)
3080 template void
3083 const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3084# endif
3085
3086# if defined(HAVE_TPETRA_INST_FLOAT)
3087 template void
3090 const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3091# endif
3092# endif
3093
3094 template void
3097 const ::LinearAlgebra::EpetraWrappers::Vector &) const;
3098
3099 template void
3100 SparseMatrix::vmult_add(MPI::Vector &, const MPI::Vector &) const;
3101
3102 template void
3104 const ::Vector<double> &) const;
3105
3106 template void
3109 const ::LinearAlgebra::distributed::Vector<double> &) const;
3110
3111# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3112# if defined(HAVE_TPETRA_INST_DOUBLE)
3113 template void
3116 const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3117# endif
3118
3119# if defined(HAVE_TPETRA_INST_FLOAT)
3120 template void
3123 const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3124# endif
3125# endif
3126
3127 template void
3130 const ::LinearAlgebra::EpetraWrappers::Vector &) const;
3131
3132 template void
3133 SparseMatrix::Tvmult_add(MPI::Vector &, const MPI::Vector &) const;
3134
3135 template void
3137 const ::Vector<double> &) const;
3138
3139 template void
3142 const ::LinearAlgebra::distributed::Vector<double> &) const;
3143
3144# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3145# if defined(HAVE_TPETRA_INST_DOUBLE)
3146 template void
3149 const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3150# endif
3151
3152# if defined(HAVE_TPETRA_INST_FLOAT)
3153 template void
3156 const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3157# endif
3158# endif
3159
3160 template void
3163 const ::LinearAlgebra::EpetraWrappers::Vector &) const;
3164} // namespace TrilinosWrappers
3165# endif // DOXYGEN
3166
3168
3169#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
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:4614
Point< 2 > first
Definition grid_out.cc:4613
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:981
@ 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:1029
Definition types.h:32
unsigned int global_dof_index
Definition types.h:81
double TrilinosScalar
Definition types.h:178