18#ifdef DEAL_II_WITH_TRILINOS
30# include <boost/container/small_vector.hpp>
32# ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT
33# include <EpetraExt_MatrixMatrix.h>
35# include <Epetra_Export.h>
36# include <Teuchos_RCP.hpp>
37# include <ml_epetra_utils.h>
38# include <ml_struct.h>
48 template <
typename VectorType>
49 typename VectorType::value_type *
55 template <
typename VectorType>
56 const typename VectorType::value_type *
62 template <
typename VectorType>
63 typename VectorType::value_type *
69 template <
typename VectorType>
70 const typename VectorType::value_type *
71 end(
const VectorType &V)
104# ifdef DEAL_II_TRILINOS_WITH_TPETRA
105 template <
typename Number>
112 template <
typename Number>
119 template <
typename Number>
127 template <
typename Number>
164 value_cache = std::make_shared<std::vector<TrilinosScalar>>(colnums);
165 colnum_cache = std::make_shared<std::vector<size_type>>(colnums);
205 : column_space_map(new Epetra_Map(0, 0,
Utilities::Trilinos::comm_self()))
207 new Epetra_FECrsMatrix(View, *column_space_map, *column_space_map, 0))
218 const unsigned int n_max_entries_per_row)
233 matrix(new Epetra_FECrsMatrix(
239 n_max_entries_per_row,
249 const std::vector<unsigned int> &n_entries_per_row)
254 , matrix(new Epetra_FECrsMatrix(
260 reinterpret_cast<
int *>(
261 const_cast<unsigned
int *>(n_entries_per_row.data())),
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,
276 n_max_entries_per_row,
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,
291 reinterpret_cast<
int *>(
292 const_cast<unsigned
int *>(
293 n_entries_per_row.data())),
302 const IndexSet &col_parallel_partitioning,
305 : column_space_map(new Epetra_Map(
306 col_parallel_partitioning.make_trilinos_map(communicator, false)))
307 , matrix(new Epetra_FECrsMatrix(
309 row_parallel_partitioning.make_trilinos_map(communicator, false),
310 n_max_entries_per_row,
319 const IndexSet &col_parallel_partitioning,
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(
326 row_parallel_partitioning.make_trilinos_map(communicator, false),
327 reinterpret_cast<
int *>(
328 const_cast<unsigned
int *>(n_entries_per_row.data())),
337 : column_space_map(new Epetra_Map(sparsity_pattern.domain_partitioner()))
339 new Epetra_FECrsMatrix(Copy,
340 sparsity_pattern.trilinos_sparsity_pattern(),
347 "The Trilinos sparsity pattern has not been compressed."));
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)
361 other.last_action = Zero;
362 other.compressed =
false;
379 bool needs_deep_copy =
384 if (!needs_deep_copy)
391 const int row_local =
matrix->RowMap().LID(
395 int n_entries, rhs_n_entries;
397 int *index_ptr, *rhs_index_ptr;
398 int ierr = rhs.
matrix->ExtractMyRowView(row_local,
405 ierr =
matrix->ExtractMyRowView(row_local,
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)
416 needs_deep_copy =
true;
420 for (
int i = 0; i < n_entries; ++i)
421 value_ptr[i] = rhs_value_ptr[i];
431 matrix = std::make_unique<Epetra_FECrsMatrix>(*rhs.
matrix);
438 std::make_unique<Epetra_CrsMatrix>(Copy, rhs.
nonlocal_matrix->Graph());
445 template <
typename SparsityPatternType>
447 reinit_matrix(
const IndexSet &row_parallel_partitioning,
448 const IndexSet &column_parallel_partitioning,
449 const SparsityPatternType &sparsity_pattern,
450 const bool exchange_data,
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)
459 nonlocal_matrix.reset();
460 nonlocal_matrix_exporter.reset();
462 column_space_map = std::make_unique<Epetra_Map>(
465 if (column_space_map->Comm().MyPID() == 0)
468 row_parallel_partitioning.
size());
470 column_parallel_partitioning.
size());
473 Epetra_Map row_space_map =
483 trilinos_sparsity.
reinit(row_parallel_partitioning,
484 column_parallel_partitioning,
488 matrix = std::make_unique<Epetra_FECrsMatrix>(
489 Copy, trilinos_sparsity.trilinos_sparsity_pattern(),
false);
499 std::vector<int> n_entries_per_row(last_row - first_row);
502 n_entries_per_row[row - first_row] = sparsity_pattern.row_length(row);
518 std::unique_ptr<Epetra_CrsGraph> graph;
519 if (row_space_map.Comm().NumProc() > 1)
520 graph = std::make_unique<Epetra_CrsGraph>(Copy,
522 n_entries_per_row.data(),
525 graph = std::make_unique<Epetra_CrsGraph>(Copy,
528 n_entries_per_row.data(),
536 std::vector<TrilinosWrappers::types::int_type> row_indices;
540 const int row_length = sparsity_pattern.row_length(row);
544 row_indices.resize(row_length, -1);
546 typename SparsityPatternType::iterator
p =
547 sparsity_pattern.begin(row);
549 p != sparsity_pattern.end(row);
551 row_indices[col] =
p->column();
553 graph->Epetra_CrsGraph::InsertGlobalIndices(row,
562 graph->FillComplete(*column_space_map, row_space_map);
563 graph->OptimizeStorage();
571 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph,
false);
583 class Epetra_CrsGraphMod :
public Epetra_CrsGraph
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)
592 SetIndicesAreGlobal()
594 this->Epetra_CrsGraph::SetIndicesAreGlobal(
true);
604 reinit_matrix(
const IndexSet &row_parallel_partitioning,
605 const IndexSet &column_parallel_partitioning,
607 const bool exchange_data,
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)
615 nonlocal_matrix.reset();
616 nonlocal_matrix_exporter.reset();
618 column_space_map = std::make_unique<Epetra_Map>(
622 row_parallel_partitioning.
size());
624 column_parallel_partitioning.
size());
626 Epetra_Map row_space_map =
631 if (relevant_rows.size() == 0)
633 relevant_rows.set_size(
635 relevant_rows.add_range(
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"));
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(
659 row_space_map.Comm());
660 return !relevant_map.SameAs(row_space_map);
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;
668 for (
const auto global_row : relevant_rows)
670 if (row_space_map.MyGID(
672 n_entries_per_row[own++] =
674 else if (sparsity_pattern.
row_length(global_row) > 0)
676 ghost_rows.push_back(global_row);
677 n_entries_per_ghost_row.push_back(
683 Epetra_Map off_processor_map(-1,
685 (ghost_rows.size() > 0) ?
686 (ghost_rows.data()) :
689 row_space_map.Comm());
691 std::unique_ptr<Epetra_CrsGraph> graph;
692 std::unique_ptr<Epetra_CrsGraphMod> nonlocal_graph;
693 if (row_space_map.Comm().NumProc() > 1)
696 std::make_unique<Epetra_CrsGraph>(Copy,
698 (n_entries_per_row.size() > 0) ?
699 (n_entries_per_row.data()) :
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());
708 std::make_unique<Epetra_CrsGraph>(Copy,
711 (n_entries_per_row.size() > 0) ?
712 (n_entries_per_row.data()) :
717 std::vector<TrilinosWrappers::types::int_type> row_indices;
719 for (
const auto global_row : relevant_rows)
721 const int row_length = sparsity_pattern.
row_length(global_row);
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);
729 if (row_space_map.MyGID(
731 graph->InsertGlobalIndices(global_row,
737 nonlocal_graph->InsertGlobalIndices(global_row,
744 if (nonlocal_graph.get() !=
nullptr)
750 nonlocal_graph->SetIndicesAreGlobal();
751 Assert(nonlocal_graph->IndicesAreGlobal() ==
true,
753 nonlocal_graph->FillComplete(*column_space_map, row_space_map);
754 nonlocal_graph->OptimizeStorage();
759 Epetra_Export exporter(nonlocal_graph->RowMap(), row_space_map);
760 int ierr = graph->Export(*nonlocal_graph, exporter, Add);
766 std::make_unique<Epetra_CrsMatrix>(Copy, *nonlocal_graph);
769 graph->FillComplete(*column_space_map, row_space_map);
770 graph->OptimizeStorage();
775 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph,
false);
781 template <
typename SparsityPatternType>
798 template <
typename SparsityPatternType>
799 inline std::enable_if_t<
800 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
802 const IndexSet &col_parallel_partitioning,
803 const SparsityPatternType &sparsity_pattern,
805 const bool exchange_data)
807 reinit_matrix(row_parallel_partitioning,
808 col_parallel_partitioning,
834 matrix = std::make_unique<Epetra_FECrsMatrix>(
839 std::make_unique<Epetra_CrsMatrix>(Copy,
853 if (
this == &sparse_matrix)
857 std::make_unique<Epetra_Map>(sparse_matrix.
trilinos_matrix().DomainMap());
860 matrix = std::make_unique<Epetra_FECrsMatrix>(
875 template <
typename number>
878 const IndexSet &row_parallel_partitioning,
879 const IndexSet &col_parallel_partitioning,
880 const ::SparseMatrix<number> &dealii_sparse_matrix,
882 const double drop_tolerance,
883 const bool copy_values,
884 const ::SparsityPattern *use_this_sparsity)
886 if (copy_values ==
false)
890 if (use_this_sparsity ==
nullptr)
891 reinit(row_parallel_partitioning,
892 col_parallel_partitioning,
893 dealii_sparse_matrix.get_sparsity_pattern(),
897 reinit(row_parallel_partitioning,
898 col_parallel_partitioning,
905 const size_type n_rows = dealii_sparse_matrix.m();
910 const ::SparsityPattern &sparsity_pattern =
911 (use_this_sparsity !=
nullptr) ?
913 dealii_sparse_matrix.get_sparsity_pattern();
915 if (
matrix.get() ==
nullptr ||
m() != n_rows ||
918 reinit(row_parallel_partitioning,
919 col_parallel_partitioning,
931 std::vector<size_type> row_indices(maximum_row_length);
932 std::vector<TrilinosScalar> values(maximum_row_length);
934 for (
size_type row = 0; row < n_rows; ++row)
936 if (row_parallel_partitioning.
is_element(row) ==
true)
939 sparsity_pattern.begin(row);
940 typename ::SparseMatrix<number>::const_iterator it =
941 dealii_sparse_matrix.begin(row);
943 if (sparsity_pattern.n_rows() == sparsity_pattern.n_cols())
947 if (std::fabs(it->value()) > drop_tolerance)
949 values[col] = it->value();
950 row_indices[col++] = it->column();
956 while (it != dealii_sparse_matrix.end(row) &&
957 select_index != sparsity_pattern.end(row))
959 while (select_index->column() < it->column() &&
960 select_index != sparsity_pattern.end(row))
962 while (it->column() < select_index->column() &&
963 it != dealii_sparse_matrix.end(row))
966 if (it == dealii_sparse_matrix.end(row))
968 if (std::fabs(it->value()) > drop_tolerance)
970 values[col] = it->value();
971 row_indices[col++] = it->column();
978 reinterpret_cast<size_type *
>(row_indices.data()),
987 template <
typename number>
990 const ::SparseMatrix<number> &dealii_sparse_matrix,
991 const double drop_tolerance,
992 const bool copy_values,
993 const ::SparsityPattern *use_this_sparsity)
997 dealii_sparse_matrix,
1008 const bool copy_values)
1010 Assert(input_matrix.Filled() ==
true,
1011 ExcMessage(
"Input CrsMatrix has not called FillComplete()!"));
1015 const Epetra_CrsGraph *graph = &input_matrix.Graph();
1020 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph,
false);
1024 if (copy_values ==
true)
1030 const size_type my_nonzeros = input_matrix.NumMyNonzeros();
1031 std::memcpy(values, in_values, my_nonzeros *
sizeof(
TrilinosScalar));
1055 "compress() can only be called with VectorOperation add, insert, or unknown"));
1062 ExcMessage(
"Operation and argument to compress() do not match"));
1088 ierr =
matrix->OptimizeStorage();
1133 matrix->ExtractMyRowView(local_row, num_entries, values, col_indices);
1138 const std::ptrdiff_t diag_index =
1139 std::find(col_indices, col_indices + num_entries, local_row) -
1143 if (diag_index != j || new_diag_value == 0)
1146 if (diag_index != num_entries)
1147 values[diag_index] = new_diag_value;
1157 for (
const auto row : rows)
1180 if (trilinos_i == -1)
1195 int nnz_present =
matrix->NumMyEntries(trilinos_i);
1204 int ierr =
matrix->ExtractMyRowView(trilinos_i,
1211 Assert(nnz_present == nnz_extracted,
1217 const std::ptrdiff_t local_col_index =
1218 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1229 if (local_col_index == nnz_present)
1234 value = values[local_col_index];
1260 if ((trilinos_i == -1) || (trilinos_j == -1))
1271 int nnz_present =
matrix->NumMyEntries(trilinos_i);
1279 int ierr =
matrix->ExtractMyRowView(trilinos_i,
1286 Assert(nnz_present == nnz_extracted,
1292 const std::ptrdiff_t local_col_index =
1293 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1303 if (local_col_index == nnz_present)
1306 value = values[local_col_index];
1352 int ierr =
matrix->NumMyRowEntries(local_row, ncols);
1356 return static_cast<unsigned int>(ncols);
1363 const std::vector<size_type> &col_indices,
1365 const bool elide_zero_values)
1367 Assert(row_indices.size() == values.m(),
1369 Assert(col_indices.size() == values.n(),
1372 for (
size_type i = 0; i < row_indices.size(); ++i)
1384 const std::vector<size_type> &col_indices,
1385 const std::vector<TrilinosScalar> &values,
1386 const bool elide_zero_values)
1388 Assert(col_indices.size() == values.size(),
1402 SparseMatrix::set<TrilinosScalar>(
const size_type row,
1406 const bool elide_zero_values)
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);
1434 if (elide_zero_values ==
false)
1439 col_value_ptr = values;
1446 col_index_ptr = local_index_array.data();
1447 col_value_ptr = local_value_array.data();
1452 const double value = values[j];
1456 local_index_array[n_columns] = col_indices[j];
1457 local_value_array[n_columns] = value;
1474 if (
matrix->RowMap().MyGID(
1477 if (
matrix->Filled() ==
false)
1479 ierr =
matrix->Epetra_CrsMatrix::InsertGlobalValues(
1480 row,
static_cast<int>(n_columns), col_value_ptr, col_index_ptr);
1489 ierr =
matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row,
1504 if (
matrix->Filled() ==
false)
1506 ierr =
matrix->InsertGlobalValues(1,
1511 Epetra_FECrsMatrix::ROW_MAJOR);
1516 ierr =
matrix->ReplaceGlobalValues(1,
1521 Epetra_FECrsMatrix::ROW_MAJOR);
1538 const bool elide_zero_values)
1540 Assert(indices.size() == values.m(),
1544 for (
size_type i = 0; i < indices.size(); ++i)
1556 const std::vector<size_type> &col_indices,
1558 const bool elide_zero_values)
1560 Assert(row_indices.size() == values.m(),
1562 Assert(col_indices.size() == values.n(),
1565 for (
size_type i = 0; i < row_indices.size(); ++i)
1577 const std::vector<size_type> &col_indices,
1578 const std::vector<TrilinosScalar> &values,
1579 const bool elide_zero_values)
1581 Assert(col_indices.size() == values.size(),
1598 const bool elide_zero_values,
1620 boost::container::small_vector<TrilinosScalar, 100> local_value_array(
1622 boost::container::small_vector<TrilinosWrappers::types::int_type, 100>
1623 local_index_array(n_cols);
1628 if (elide_zero_values ==
false)
1633 col_value_ptr = values;
1644 col_index_ptr = local_index_array.data();
1645 col_value_ptr = local_value_array.data();
1650 const double value = values[j];
1655 local_index_array[n_columns] = col_indices[j];
1656 local_value_array[n_columns] = value;
1672 if (
matrix->RowMap().MyGID(
1675 ierr =
matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row,
1689 ExcMessage(
"Attempted to write into off-processor matrix row "
1690 "that has not be specified as being writable upon "
1708 ierr =
matrix->SumIntoGlobalValues(1,
1713 Epetra_FECrsMatrix::ROW_MAJOR);
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;
1725 std::cout << col_index_ptr[i] <<
" ";
1726 std::cout << std::endl << std::endl;
1727 std::cout <<
"Matrix row "
1728 << (
matrix->RowMap().MyGID(
1733 <<
" has the following indices:" << std::endl;
1734 std::vector<TrilinosWrappers::types::int_type> indices;
1735 const Epetra_CrsGraph *graph =
1742 indices.resize(graph->NumGlobalIndices(row));
1744 graph->ExtractGlobalRowCopy(row,
1751 std::cout << indices[i] <<
" ";
1752 std::cout << std::endl << std::endl;
1773 const int ierr =
matrix->PutScalar(0.0);
1795 ExcMessage(
"Can only add matrices with same distribution of rows"));
1797 ExcMessage(
"Addition of matrices only allowed if matrices are "
1798 "filled, i.e., compress() has been called"));
1800 const bool same_col_map =
matrix->ColMap().SameAs(rhs.
matrix->ColMap());
1804 const int row_local =
matrix->RowMap().LID(
1812 int n_entries, rhs_n_entries;
1814 int *index_ptr, *rhs_index_ptr;
1815 int ierr = rhs.
matrix->ExtractMyRowView(row_local,
1823 matrix->ExtractMyRowView(row_local, n_entries, value_ptr, index_ptr);
1825 bool expensive_checks = (n_entries != rhs_n_entries || !same_col_map);
1826 if (!expensive_checks)
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;
1841 if (expensive_checks)
1843 for (
int i = 0; i < rhs_n_entries; ++i)
1845 if (rhs_value_ptr[i] == 0.)
1849 int local_col =
matrix->ColMap().LID(rhs_global_col);
1851 index_ptr + n_entries,
1853 Assert(local_index != index_ptr + n_entries &&
1854 *local_index == local_col,
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];
1878 if (!
matrix->UseTranspose())
1880 ierr =
matrix->SetUseTranspose(
true);
1885 ierr =
matrix->SetUseTranspose(
false);
1895 const int ierr =
matrix->Scale(a);
1911 const int ierr =
matrix->Scale(factor);
1924 return matrix->NormOne();
1933 return matrix->NormInf();
1942 return matrix->NormFrobenius();
1949 namespace SparseMatrixImplementation
1951 template <
typename VectorType>
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?"));
1983 template <
typename VectorType>
1985 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2001 Epetra_MultiVector tril_dst(
2003 Epetra_MultiVector tril_src(View,
2010 const int ierr =
matrix->Multiply(
false, tril_src, tril_dst);
2017 template <
typename VectorType>
2019 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2027 template <
typename VectorType>
2029 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2043 Epetra_MultiVector tril_dst(
2045 Epetra_MultiVector tril_src(View,
2051 const int ierr =
matrix->Multiply(
true, tril_src, tril_dst);
2058 template <
typename VectorType>
2060 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2068 template <
typename VectorType>
2076 VectorType tmp_vector;
2077 tmp_vector.reinit(dst,
true);
2078 vmult(tmp_vector, src);
2084 template <
typename VectorType>
2092 VectorType tmp_vector;
2093 tmp_vector.reinit(dst,
true);
2107 temp_vector.
reinit(v,
true);
2109 vmult(temp_vector, v);
2110 return temp_vector * v;
2124 temp_vector.
reinit(v,
true);
2126 vmult(temp_vector, v);
2127 return u * temp_vector;
2141 const bool transpose_left)
2143 const bool use_vector = (V.
size() == inputright.
m() ? true :
false);
2144 if (transpose_left ==
false)
2146 Assert(inputleft.
n() == inputright.
m(),
2150 ExcMessage(
"Parallel partitioning of A and B does not fit."));
2154 Assert(inputleft.
m() == inputright.
m(),
2158 ExcMessage(
"Parallel partitioning of A and B does not fit."));
2169 Teuchos::RCP<Epetra_CrsMatrix> mod_B;
2170 if (use_vector ==
false)
2172 mod_B = Teuchos::rcp(
const_cast<Epetra_CrsMatrix *
>(
2178 mod_B = Teuchos::rcp(
2184 ExcMessage(
"Parallel distribution of matrix B and vector V "
2185 "does not match."));
2188 for (
int i = 0; i < local_N; ++i)
2191 double *new_data, *B_data;
2192 mod_B->ExtractMyRowView(i, N_entries, new_data);
2198 new_data[j] = value * B_data[j];
2209# ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT
2214 const_cast<Epetra_CrsMatrix &
>(
2215 tmp_result.trilinos_matrix()));
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."));
2224 result.
reinit(tmp_result.trilinos_matrix());
2262 const bool print_detailed_trilinos_information)
const
2264 if (print_detailed_trilinos_information ==
true)
2272 for (
int i = 0; i <
matrix->NumMyRows(); ++i)
2275 matrix->ExtractMyRowView(i, num_entries, values, indices);
2283 <<
") " << values[j] << std::endl;
2296 sizeof(*this) +
sizeof(*matrix) +
sizeof(*
matrix->Graph().DataPtr());
2299 matrix->NumMyNonzeros() +
2308 const Epetra_MpiComm *mpi_comm =
2309 dynamic_cast<const Epetra_MpiComm *
>(&
matrix->RangeMap().Comm());
2311 return mpi_comm->Comm();
2320 namespace LinearOperatorImplementation
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()))
2330 ExcMessage(
"Uninitialized TrilinosPayload::vmult called "
2331 "(Default constructor)"));
2336 ExcMessage(
"Uninitialized TrilinosPayload::Tvmult called "
2337 "(Default constructor)"));
2342 ExcMessage(
"Uninitialized TrilinosPayload::inv_vmult called "
2343 "(Default constructor)"));
2348 ExcMessage(
"Uninitialized TrilinosPayload::inv_Tvmult called "
2349 "(Default constructor)"));
2359 matrix.trilinos_matrix()),
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())
2374 matrix.trilinos_matrix()),
2376 payload_exemplar.UseTranspose(),
2377 payload_exemplar.get_mpi_communicator(),
2378 payload_exemplar.locally_owned_domain_indices(),
2379 payload_exemplar.locally_owned_range_indices())
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())
2401 preconditioner.trilinos_operator(),
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())
2416 payload_exemplar.UseTranspose(),
2417 payload_exemplar.get_mpi_communicator(),
2418 payload_exemplar.locally_owned_domain_indices(),
2419 payload_exemplar.locally_owned_range_indices())
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)
2441 : use_transpose(false)
2444 communicator(first_op.communicator)
2445 , domain_map(second_op.domain_map)
2446 , range_map(first_op.range_map)
2457 tril_dst = tril_src;
2461 tril_dst = tril_src;
2465 tril_dst = tril_src;
2469 tril_dst = tril_src;
2483 const int ierr = tril_dst.PutScalar(0.0);
2489 const int ierr = tril_dst.PutScalar(0.0);
2496 ExcMessage(
"Cannot compute inverse of null operator"));
2498 const int ierr = tril_dst.PutScalar(0.0);
2505 ExcMessage(
"Cannot compute inverse of null operator"));
2507 const int ierr = tril_dst.PutScalar(0.0);
2607 return "TrilinosPayload";
2665 "Operators are set to work on incompatible IndexSets."));
2669 "Operators are set to work on incompatible IndexSets."));
2674 return_op.
vmult = [first_op, second_op](Range &tril_dst,
2675 const Domain &tril_src) {
2683 i->reinit(
IndexSet(first_op_init_map),
2688 const size_type i_local_size = i->end() - i->begin();
2692 (void)second_op_init_map;
2693 Intermediate tril_int(View,
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);
2707 return_op.
Tvmult = [first_op, second_op](Domain &tril_dst,
2708 const Range &tril_src) {
2723 i->reinit(
IndexSet(first_op_init_map),
2728 const size_type i_local_size = i->end() - i->begin();
2732 (void)second_op_init_map;
2733 Intermediate tril_int(View,
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);
2751 return_op.
inv_vmult = [first_op, second_op](Domain &tril_dst,
2752 const Range &tril_src) {
2760 i->reinit(
IndexSet(first_op_init_map),
2765 const size_type i_local_size = i->end() - i->begin();
2769 (void)second_op_init_map;
2770 Intermediate tril_int(View,
2780 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2784 return_op.
inv_Tvmult = [first_op, second_op](Range &tril_dst,
2785 const Domain &tril_src) {
2800 i->reinit(
IndexSet(first_op_init_map),
2805 const size_type i_local_size = i->end() - i->begin();
2809 (void)second_op_init_map;
2810 Intermediate tril_int(View,
2820 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2845 "Operators are set to work on incompatible IndexSets."));
2850 return_op.
vmult = [first_op, second_op](Range &tril_dst,
2851 const Domain &tril_src) {
2859 i->reinit(
IndexSet(first_op_init_map),
2864 const size_type i_local_size = i->end() - i->begin();
2868 (void)second_op_init_map;
2869 Intermediate tril_int(View,
2877 second_op.
Apply(tril_src, tril_int);
2878 first_op.
Apply(tril_int, tril_dst);
2881 return_op.
Tvmult = [first_op, second_op](Domain &tril_dst,
2882 const Range &tril_src) {
2897 i->reinit(
IndexSet(first_op_init_map),
2902 const size_type i_local_size = i->end() - i->begin();
2906 (void)second_op_init_map;
2907 Intermediate tril_int(View,
2914 first_op.
Apply(tril_src, tril_int);
2915 second_op.
Apply(tril_int, tril_dst);
2922 return_op.
inv_vmult = [first_op, second_op](Domain &tril_dst,
2923 const Range &tril_src) {
2931 i->reinit(
IndexSet(first_op_init_map),
2936 const size_type i_local_size = i->end() - i->begin();
2940 (void)second_op_init_map;
2941 Intermediate tril_int(View,
2953 return_op.
inv_Tvmult = [first_op, second_op](Range &tril_dst,
2954 const Domain &tril_src) {
2969 i->reinit(
IndexSet(first_op_init_map),
2974 const size_type i_local_size = i->end() - i->begin();
2978 (void)second_op_init_map;
2979 Intermediate tril_int(View,
3007# include "trilinos_sparse_matrix.inst"
3022 const ::SparsityPattern &,
3038 const ::Vector<double> &)
const;
3043 const ::LinearAlgebra::distributed::Vector<double> &)
const;
3045# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3046# if defined(HAVE_TPETRA_INST_DOUBLE)
3050 const ::LinearAlgebra::TpetraWrappers::Vector<double> &)
const;
3053# if defined(HAVE_TPETRA_INST_FLOAT)
3057 const ::LinearAlgebra::TpetraWrappers::Vector<float> &)
const;
3064 const ::LinearAlgebra::EpetraWrappers::Vector &)
const;
3071 const ::Vector<double> &)
const;
3076 const ::LinearAlgebra::distributed::Vector<double> &)
const;
3078# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3079# if defined(HAVE_TPETRA_INST_DOUBLE)
3083 const ::LinearAlgebra::TpetraWrappers::Vector<double> &)
const;
3086# if defined(HAVE_TPETRA_INST_FLOAT)
3090 const ::LinearAlgebra::TpetraWrappers::Vector<float> &)
const;
3097 const ::LinearAlgebra::EpetraWrappers::Vector &)
const;
3104 const ::Vector<double> &)
const;
3109 const ::LinearAlgebra::distributed::Vector<double> &)
const;
3111# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3112# if defined(HAVE_TPETRA_INST_DOUBLE)
3116 const ::LinearAlgebra::TpetraWrappers::Vector<double> &)
const;
3119# if defined(HAVE_TPETRA_INST_FLOAT)
3123 const ::LinearAlgebra::TpetraWrappers::Vector<float> &)
const;
3130 const ::LinearAlgebra::EpetraWrappers::Vector &)
const;
3137 const ::Vector<double> &)
const;
3142 const ::LinearAlgebra::distributed::Vector<double> &)
const;
3144# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3145# if defined(HAVE_TPETRA_INST_DOUBLE)
3149 const ::LinearAlgebra::TpetraWrappers::Vector<double> &)
const;
3152# if defined(HAVE_TPETRA_INST_FLOAT)
3156 const ::LinearAlgebra::TpetraWrappers::Vector<float> &)
const;
3163 const ::LinearAlgebra::EpetraWrappers::Vector &)
const;
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
bool is_element(const size_type index) const
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
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
MPI_Comm get_mpi_communicator() const
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
TrilinosScalar l1_norm() const
void compress(VectorOperation::values operation)
std::unique_ptr< Epetra_FECrsMatrix > matrix
TrilinosScalar linfty_norm() const
const Epetra_CrsMatrix & trilinos_matrix() const
size_type memory_consumption() const
TrilinosScalar matrix_norm_square(const MPI::Vector &v) const
Epetra_CombineMode last_action
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
TrilinosScalar frobenius_norm() 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_Map & domain_partitioner() const
const Epetra_FECrsGraph & trilinos_sparsity_pattern() const
Epetra_MpiComm communicator
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
virtual int SetUseTranspose(bool UseTranspose) override
virtual bool UseTranspose() const override
virtual const Epetra_Map & OperatorDomainMap() const override
IndexSet locally_owned_range_indices() const
TrilinosPayload transpose_payload() const
Epetra_MultiVector VectorType
MPI_Comm get_mpi_communicator() const
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
std::function< void(VectorType &, const VectorType &)> Tvmult
virtual int Apply(const VectorType &X, VectorType &Y) const override
std::function< void(VectorType &, const VectorType &)> vmult
virtual const Epetra_Map & OperatorRangeMap() const override
std::function< void(VectorType &, const VectorType &)> inv_vmult
virtual const char * Label() const override
IndexSet locally_owned_domain_indices() const
virtual const Epetra_Comm & Comm() const override
virtual bool HasNormInf() const override
TrilinosPayload identity_payload() const
virtual double NormInf() const override
TrilinosPayload null_payload() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
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)
@ 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)
unsigned int global_dof_index