18#ifdef DEAL_II_WITH_TRILINOS
30# include <boost/container/small_vector.hpp>
33# ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT
34# include <EpetraExt_MatrixMatrix.h>
36# include <Epetra_Export.h>
37# include <Teuchos_RCP.hpp>
38# include <ml_epetra_utils.h>
39# include <ml_struct.h>
50 template <
typename VectorType>
51 typename VectorType::value_type *
57 template <
typename VectorType>
58 const typename VectorType::value_type *
64 template <
typename VectorType>
65 typename VectorType::value_type *
71 template <
typename VectorType>
72 const typename VectorType::value_type *
73 end(
const VectorType &V)
82 return V.trilinos_vector()[0];
89 return V.trilinos_vector()[0];
96 return V.trilinos_vector()[0] + V.trilinos_vector().MyLength();
103 return V.trilinos_vector()[0] + V.trilinos_vector().MyLength();
106# ifdef DEAL_II_TRILINOS_WITH_TPETRA
107 template <
typename Number,
typename MemorySpace>
111 return V.trilinos_vector().getDataNonConst().get();
114 template <
typename Number,
typename MemorySpace>
118 return V.trilinos_vector().getData().get();
121 template <
typename Number,
typename MemorySpace>
125 return V.trilinos_vector().getDataNonConst().get() +
126 V.trilinos_vector().getLocalLength();
129 template <
typename Number,
typename MemorySpace>
133 return V.trilinos_vector().getData().get() +
134 V.trilinos_vector().getLocalLength();
166 value_cache = std::make_shared<std::vector<TrilinosScalar>>(colnums);
167 colnum_cache = std::make_shared<std::vector<size_type>>(colnums);
207 : column_space_map(new Epetra_Map(0, 0,
Utilities::Trilinos::comm_self()))
209 new Epetra_FECrsMatrix(View, *column_space_map, *column_space_map, 0))
220 const unsigned int n_max_entries_per_row)
235 matrix(new Epetra_FECrsMatrix(
241 n_max_entries_per_row,
251 const std::vector<unsigned int> &n_entries_per_row)
256 , matrix(new Epetra_FECrsMatrix(
262 reinterpret_cast<
int *>(
263 const_cast<unsigned
int *>(n_entries_per_row.
data())),
273 const unsigned int n_max_entries_per_row)
274 : column_space_map(new Epetra_Map(
275 parallel_partitioning.make_trilinos_map(communicator, false)))
276 , matrix(new Epetra_FECrsMatrix(Copy,
278 n_max_entries_per_row,
288 const std::vector<unsigned int> &n_entries_per_row)
289 : column_space_map(new Epetra_Map(
290 parallel_partitioning.make_trilinos_map(communicator, false)))
291 , matrix(new Epetra_FECrsMatrix(Copy,
293 reinterpret_cast<
int *>(
294 const_cast<unsigned
int *>(
295 n_entries_per_row.
data())),
304 const IndexSet &col_parallel_partitioning,
307 : column_space_map(new Epetra_Map(
308 col_parallel_partitioning.make_trilinos_map(communicator, false)))
309 , matrix(new Epetra_FECrsMatrix(
311 row_parallel_partitioning.make_trilinos_map(communicator, false),
312 n_max_entries_per_row,
321 const IndexSet &col_parallel_partitioning,
323 const std::vector<unsigned int> &n_entries_per_row)
324 : column_space_map(new Epetra_Map(
325 col_parallel_partitioning.make_trilinos_map(communicator, false)))
326 , matrix(new Epetra_FECrsMatrix(
328 row_parallel_partitioning.make_trilinos_map(communicator, false),
329 reinterpret_cast<
int *>(
330 const_cast<unsigned
int *>(n_entries_per_row.
data())),
339 : column_space_map(new Epetra_Map(sparsity_pattern.domain_partitioner()))
341 new Epetra_FECrsMatrix(Copy,
342 sparsity_pattern.trilinos_sparsity_pattern(),
349 "The Trilinos sparsity pattern has not been compressed."));
356 : column_space_map(std::move(other.column_space_map))
357 , matrix(std::move(other.matrix))
358 , nonlocal_matrix(std::move(other.nonlocal_matrix))
359 , nonlocal_matrix_exporter(std::move(other.nonlocal_matrix_exporter))
360 , last_action(other.last_action)
361 , compressed(other.compressed)
363 other.last_action = Zero;
364 other.compressed =
false;
381 bool needs_deep_copy =
386 if (!needs_deep_copy)
393 const int row_local =
matrix->RowMap().LID(
397 int n_entries, rhs_n_entries;
399 int *index_ptr, *rhs_index_ptr;
400 int ierr = rhs.
matrix->ExtractMyRowView(row_local,
406 ierr =
matrix->ExtractMyRowView(row_local,
412 if (n_entries != rhs_n_entries ||
413 std::memcmp(
static_cast<void *
>(index_ptr),
414 static_cast<void *
>(rhs_index_ptr),
415 sizeof(
int) * n_entries) != 0)
417 needs_deep_copy =
true;
421 for (
int i = 0; i < n_entries; ++i)
422 value_ptr[i] = rhs_value_ptr[i];
432 matrix = std::make_unique<Epetra_FECrsMatrix>(*rhs.
matrix);
439 std::make_unique<Epetra_CrsMatrix>(Copy, rhs.
nonlocal_matrix->Graph());
446 template <
typename SparsityPatternType>
448 reinit_matrix(
const IndexSet &row_parallel_partitioning,
449 const IndexSet &column_parallel_partitioning,
450 const SparsityPatternType &sparsity_pattern,
451 const bool exchange_data,
453 std::unique_ptr<Epetra_Map> &column_space_map,
454 std::unique_ptr<Epetra_FECrsMatrix> &matrix,
455 std::unique_ptr<Epetra_CrsMatrix> &nonlocal_matrix,
456 std::unique_ptr<Epetra_Export> &nonlocal_matrix_exporter)
460 nonlocal_matrix.reset();
461 nonlocal_matrix_exporter.reset();
463 column_space_map = std::make_unique<Epetra_Map>(
466 if (column_space_map->Comm().MyPID() == 0)
469 row_parallel_partitioning.
size());
471 column_parallel_partitioning.
size());
474 Epetra_Map row_space_map =
484 trilinos_sparsity.
reinit(row_parallel_partitioning,
485 column_parallel_partitioning,
489 matrix = std::make_unique<Epetra_FECrsMatrix>(
490 Copy, trilinos_sparsity.trilinos_sparsity_pattern(),
false);
500 std::vector<int> n_entries_per_row(last_row - first_row);
503 n_entries_per_row[row - first_row] = sparsity_pattern.row_length(row);
519 std::unique_ptr<Epetra_CrsGraph> graph;
520 if (row_space_map.Comm().NumProc() > 1)
521 graph = std::make_unique<Epetra_CrsGraph>(Copy,
523 n_entries_per_row.data(),
526 graph = std::make_unique<Epetra_CrsGraph>(Copy,
529 n_entries_per_row.data(),
537 std::vector<TrilinosWrappers::types::int_type> row_indices;
541 const int row_length = sparsity_pattern.row_length(row);
545 row_indices.resize(row_length, -1);
547 typename SparsityPatternType::iterator p =
548 sparsity_pattern.begin(row);
550 p != sparsity_pattern.end(row);
552 row_indices[col] = p->column();
554 graph->Epetra_CrsGraph::InsertGlobalIndices(row,
563 graph->FillComplete(*column_space_map, row_space_map);
564 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);
765 std::make_unique<Epetra_CrsMatrix>(Copy, *nonlocal_graph);
768 graph->FillComplete(*column_space_map, row_space_map);
769 graph->OptimizeStorage();
774 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph,
false);
780 template <
typename SparsityPatternType>
797 template <
typename SparsityPatternType>
798 inline std::enable_if_t<
799 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
801 const IndexSet &col_parallel_partitioning,
802 const SparsityPatternType &sparsity_pattern,
804 const bool exchange_data)
806 reinit_matrix(row_parallel_partitioning,
807 col_parallel_partitioning,
833 matrix = std::make_unique<Epetra_FECrsMatrix>(
838 std::make_unique<Epetra_CrsMatrix>(Copy,
852 if (
this == &sparse_matrix)
856 std::make_unique<Epetra_Map>(sparse_matrix.
trilinos_matrix().DomainMap());
859 matrix = std::make_unique<Epetra_FECrsMatrix>(
874 template <
typename number>
877 const IndexSet &row_parallel_partitioning,
878 const IndexSet &col_parallel_partitioning,
879 const ::SparseMatrix<number> &dealii_sparse_matrix,
881 const double drop_tolerance,
882 const bool copy_values,
883 const ::SparsityPattern *use_this_sparsity)
885 if (copy_values ==
false)
889 if (use_this_sparsity ==
nullptr)
890 reinit(row_parallel_partitioning,
891 col_parallel_partitioning,
892 dealii_sparse_matrix.get_sparsity_pattern(),
896 reinit(row_parallel_partitioning,
897 col_parallel_partitioning,
904 const size_type n_rows = dealii_sparse_matrix.m();
909 const ::SparsityPattern &sparsity_pattern =
910 (use_this_sparsity !=
nullptr) ?
912 dealii_sparse_matrix.get_sparsity_pattern();
914 if (
matrix.get() ==
nullptr ||
m() != n_rows ||
917 reinit(row_parallel_partitioning,
918 col_parallel_partitioning,
930 std::vector<size_type> row_indices(maximum_row_length);
931 std::vector<TrilinosScalar> values(maximum_row_length);
933 for (
size_type row = 0; row < n_rows; ++row)
935 if (row_parallel_partitioning.
is_element(row) ==
true)
938 sparsity_pattern.begin(row);
939 typename ::SparseMatrix<number>::const_iterator it =
940 dealii_sparse_matrix.begin(row);
942 if (sparsity_pattern.n_rows() == sparsity_pattern.n_cols())
946 if (std::fabs(it->value()) > drop_tolerance)
948 values[col] = it->value();
949 row_indices[col++] = it->column();
955 while (it != dealii_sparse_matrix.end(row) &&
956 select_index != sparsity_pattern.end(row))
958 while (select_index->column() < it->column() &&
959 select_index != sparsity_pattern.end(row))
961 while (it->column() < select_index->column() &&
962 it != dealii_sparse_matrix.end(row))
965 if (it == dealii_sparse_matrix.end(row))
967 if (std::fabs(it->value()) > drop_tolerance)
969 values[col] = it->value();
970 row_indices[col++] = it->column();
977 reinterpret_cast<size_type *
>(row_indices.data()),
986 template <
typename number>
989 const ::SparseMatrix<number> &dealii_sparse_matrix,
990 const double drop_tolerance,
991 const bool copy_values,
992 const ::SparsityPattern *use_this_sparsity)
996 dealii_sparse_matrix,
1007 const bool copy_values)
1009 Assert(input_matrix.Filled() ==
true,
1010 ExcMessage(
"Input CrsMatrix has not called FillComplete()!"));
1014 const Epetra_CrsGraph *graph = &input_matrix.Graph();
1019 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph,
false);
1023 if (copy_values ==
true)
1029 const size_type my_nonzeros = input_matrix.NumMyNonzeros();
1030 std::memcpy(values, in_values, my_nonzeros *
sizeof(
TrilinosScalar));
1054 "compress() can only be called with VectorOperation add, insert, or unknown"));
1061 ExcMessage(
"Operation and argument to compress() do not match"));
1087 ierr =
matrix->OptimizeStorage();
1132 matrix->ExtractMyRowView(local_row, num_entries, values, col_indices);
1136 const std::ptrdiff_t diag_index =
1137 std::find(col_indices, col_indices + num_entries, local_row) -
1141 if (diag_index != j || new_diag_value == 0)
1144 if (diag_index != num_entries)
1145 values[diag_index] = new_diag_value;
1155 for (
const auto row : rows)
1178 if (trilinos_i == -1)
1193 int nnz_present =
matrix->NumMyEntries(trilinos_i);
1202 int ierr =
matrix->ExtractMyRowView(trilinos_i,
1208 Assert(nnz_present == nnz_extracted,
1214 const std::ptrdiff_t local_col_index =
1215 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1226 if (local_col_index == nnz_present)
1231 value = values[local_col_index];
1257 if ((trilinos_i == -1) || (trilinos_j == -1))
1268 int nnz_present =
matrix->NumMyEntries(trilinos_i);
1276 int ierr =
matrix->ExtractMyRowView(trilinos_i,
1282 Assert(nnz_present == nnz_extracted,
1288 const std::ptrdiff_t local_col_index =
1289 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1299 if (local_col_index == nnz_present)
1302 value = values[local_col_index];
1351 int ierr =
matrix->NumMyRowEntries(local_row, ncols);
1355 return static_cast<unsigned int>(ncols);
1362 const std::vector<size_type> &col_indices,
1364 const bool elide_zero_values)
1366 Assert(row_indices.size() == values.m(),
1368 Assert(col_indices.size() == values.n(),
1371 for (
size_type i = 0; i < row_indices.size(); ++i)
1383 const std::vector<size_type> &col_indices,
1384 const std::vector<TrilinosScalar> &values,
1385 const bool elide_zero_values)
1387 Assert(col_indices.size() == values.size(),
1401 SparseMatrix::set<TrilinosScalar>(
const size_type row,
1405 const bool elide_zero_values)
1425 boost::container::small_vector<TrilinosScalar, 200> local_value_array(
1426 elide_zero_values ? n_cols : 0);
1427 boost::container::small_vector<TrilinosWrappers::types::int_type, 200>
1428 local_index_array(elide_zero_values ? n_cols : 0);
1433 if (elide_zero_values ==
false)
1438 col_value_ptr = values;
1445 col_index_ptr = local_index_array.data();
1446 col_value_ptr = local_value_array.data();
1451 const double value = values[j];
1455 local_index_array[n_columns] = col_indices[j];
1456 local_value_array[n_columns] = value;
1473 if (
matrix->RowMap().MyGID(
1476 if (
matrix->Filled() ==
false)
1478 ierr =
matrix->Epetra_CrsMatrix::InsertGlobalValues(
1479 row,
static_cast<int>(n_columns), col_value_ptr, col_index_ptr);
1488 ierr =
matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row,
1503 if (
matrix->Filled() ==
false)
1505 ierr =
matrix->InsertGlobalValues(1,
1510 Epetra_FECrsMatrix::ROW_MAJOR);
1515 ierr =
matrix->ReplaceGlobalValues(1,
1520 Epetra_FECrsMatrix::ROW_MAJOR);
1537 const bool elide_zero_values)
1539 Assert(indices.size() == values.m(),
1543 for (
size_type i = 0; i < indices.size(); ++i)
1555 const std::vector<size_type> &col_indices,
1557 const bool elide_zero_values)
1559 Assert(row_indices.size() == values.m(),
1561 Assert(col_indices.size() == values.n(),
1564 for (
size_type i = 0; i < row_indices.size(); ++i)
1576 const std::vector<size_type> &col_indices,
1577 const std::vector<TrilinosScalar> &values,
1578 const bool elide_zero_values)
1580 Assert(col_indices.size() == values.size(),
1597 const bool elide_zero_values,
1622 boost::container::small_vector<TrilinosScalar, 100> local_value_array(
1624 boost::container::small_vector<TrilinosWrappers::types::int_type, 100>
1625 local_index_array(n_cols);
1630 if (elide_zero_values ==
false)
1635 col_value_ptr = values;
1647 col_index_ptr = local_index_array.data();
1648 col_value_ptr = local_value_array.data();
1653 const double value = values[j];
1658 local_index_array[n_columns] = col_indices[j];
1659 local_value_array[n_columns] = value;
1675 if (
matrix->RowMap().MyGID(
1678 ierr =
matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row,
1692 ExcMessage(
"Attempted to write into off-processor matrix row "
1693 "that has not be specified as being writable upon "
1711 ierr =
matrix->SumIntoGlobalValues(1,
1716 Epetra_FECrsMatrix::ROW_MAJOR);
1724 std::cout <<
"------------------------------------------"
1726 std::cout <<
"Got error " << ierr <<
" in row " << row
1727 <<
" of proc " <<
matrix->RowMap().Comm().MyPID()
1728 <<
" when trying to add the columns:" << std::endl;
1730 std::cout << col_index_ptr[i] <<
" ";
1731 std::cout << std::endl << std::endl;
1732 std::cout <<
"Matrix row "
1733 << (
matrix->RowMap().MyGID(
1738 <<
" has the following indices:" << std::endl;
1739 std::vector<TrilinosWrappers::types::int_type> indices;
1740 const Epetra_CrsGraph *graph =
1748 indices.resize(graph->NumGlobalIndices(row));
1750 graph->ExtractGlobalRowCopy(row,
1757 std::cout << indices[i] <<
" ";
1758 std::cout << std::endl << std::endl;
1779 const int ierr =
matrix->PutScalar(0.0);
1801 ExcMessage(
"Can only add matrices with same distribution of rows"));
1803 ExcMessage(
"Addition of matrices only allowed if matrices are "
1804 "filled, i.e., compress() has been called"));
1806 const bool same_col_map =
matrix->ColMap().SameAs(rhs.
matrix->ColMap());
1810 const int row_local =
matrix->RowMap().LID(
1818 int n_entries, rhs_n_entries;
1820 int *index_ptr, *rhs_index_ptr;
1821 int ierr = rhs.
matrix->ExtractMyRowView(row_local,
1828 matrix->ExtractMyRowView(row_local, n_entries, value_ptr, index_ptr);
1830 bool expensive_checks = (n_entries != rhs_n_entries || !same_col_map);
1831 if (!expensive_checks)
1835 expensive_checks = std::memcmp(
static_cast<void *
>(index_ptr),
1836 static_cast<void *
>(rhs_index_ptr),
1837 sizeof(
int) * n_entries) != 0;
1838 if (!expensive_checks)
1839 for (
int i = 0; i < n_entries; ++i)
1840 value_ptr[i] += rhs_value_ptr[i] * factor;
1846 if (expensive_checks)
1848 for (
int i = 0; i < rhs_n_entries; ++i)
1850 if (rhs_value_ptr[i] == 0.)
1854 int local_col =
matrix->ColMap().LID(rhs_global_col);
1856 index_ptr + n_entries,
1858 Assert(local_index != index_ptr + n_entries &&
1859 *local_index == local_col,
1861 "Adding the entries from the other matrix "
1862 "failed, because the sparsity pattern "
1863 "of that matrix includes more elements than the "
1864 "calling matrix, which is not allowed."));
1865 value_ptr[local_index - index_ptr] += factor * rhs_value_ptr[i];
1883 if (!
matrix->UseTranspose())
1885 ierr =
matrix->SetUseTranspose(
true);
1890 ierr =
matrix->SetUseTranspose(
false);
1900 const int ierr =
matrix->Scale(a);
1915 const int ierr =
matrix->Scale(factor);
1927 return matrix->NormOne();
1936 return matrix->NormInf();
1945 return matrix->NormFrobenius();
1952 namespace SparseMatrixImplementation
1954 template <
typename VectorType>
1967 ExcMessage(
"The column partitioning of a matrix does not match "
1968 "the partitioning of a vector you are trying to "
1969 "multiply it with. Are you multiplying the "
1970 "matrix with a vector that has ghost elements?"));
1972 ExcMessage(
"The row partitioning of a matrix does not match "
1973 "the partitioning of a vector you are trying to "
1974 "put the result of a matrix-vector product in. "
1975 "Are you trying to put the product of the "
1976 "matrix with a vector into a vector that has "
1977 "ghost elements?"));
1983 template <
typename VectorType>
1985 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1999 Epetra_MultiVector tril_dst(
2001 Epetra_MultiVector tril_src(View,
2008 const int ierr =
matrix->Multiply(
false, tril_src, tril_dst);
2014 template <
typename VectorType>
2016 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2024 template <
typename VectorType>
2026 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2040 Epetra_MultiVector tril_dst(
2042 Epetra_MultiVector tril_src(View,
2048 const int ierr =
matrix->Multiply(
true, tril_src, tril_dst);
2054 template <
typename VectorType>
2056 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
2064 template <
typename VectorType>
2072 VectorType tmp_vector;
2073 tmp_vector.reinit(dst,
true);
2074 vmult(tmp_vector, src);
2080 template <
typename VectorType>
2088 VectorType tmp_vector;
2089 tmp_vector.reinit(dst,
true);
2103 temp_vector.
reinit(v,
true);
2105 vmult(temp_vector, v);
2106 return temp_vector * v;
2120 temp_vector.
reinit(v,
true);
2122 vmult(temp_vector, v);
2123 return u * temp_vector;
2135 const bool transpose_left)
2137 const bool use_vector = (V.size() == inputright.
m() ? true :
false);
2138 if (transpose_left ==
false)
2140 Assert(inputleft.
n() == inputright.
m(),
2144 ExcMessage(
"Parallel partitioning of A and B does not fit."));
2148 Assert(inputleft.
m() == inputright.
m(),
2152 ExcMessage(
"Parallel partitioning of A and B does not fit."));
2163 Teuchos::RCP<Epetra_CrsMatrix> mod_B;
2166 mod_B = Teuchos::rcp(
const_cast<Epetra_CrsMatrix *
>(
2172 mod_B = Teuchos::rcp(
2178 ExcMessage(
"Parallel distribution of matrix B and vector V "
2179 "does not match."));
2182 for (
int i = 0; i < local_N; ++i)
2185 double *new_data, *B_data;
2186 mod_B->ExtractMyRowView(i, N_entries, new_data);
2190 double value = V.trilinos_vector()[0][i];
2192 new_data[j] = value * B_data[j];
2203# ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT
2208 const_cast<Epetra_CrsMatrix &
>(
2209 tmp_result.trilinos_matrix()));
2212 ExcMessage(
"This function requires that the Trilinos "
2213 "installation found while running the deal.II "
2214 "CMake scripts contains the optional Trilinos "
2215 "package 'EpetraExt'. However, this optional "
2216 "part of Trilinos was not found."));
2218 result.
reinit(tmp_result.trilinos_matrix());
2256 const bool print_detailed_trilinos_information)
const
2258 if (print_detailed_trilinos_information ==
true)
2266 for (
int i = 0; i <
matrix->NumMyRows(); ++i)
2269 matrix->ExtractMyRowView(i, num_entries, values, indices);
2276 <<
") " << values[j] << std::endl;
2289 sizeof(*this) +
sizeof(*matrix) +
sizeof(*
matrix->Graph().DataPtr());
2292 matrix->NumMyNonzeros() +
2301 const Epetra_MpiComm *mpi_comm =
2302 dynamic_cast<const Epetra_MpiComm *
>(&
matrix->RangeMap().Comm());
2304 return mpi_comm->Comm();
2313 namespace LinearOperatorImplementation
2318 : use_transpose(false)
2319 , communicator(MPI_COMM_SELF)
2320 , domain_map(
IndexSet().make_trilinos_map(communicator.Comm()))
2321 , range_map(
IndexSet().make_trilinos_map(communicator.Comm()))
2325 ExcMessage(
"Uninitialized TrilinosPayload::vmult called "
2326 "(Default constructor)"));
2331 ExcMessage(
"Uninitialized TrilinosPayload::Tvmult called "
2332 "(Default constructor)"));
2337 ExcMessage(
"Uninitialized TrilinosPayload::inv_vmult called "
2338 "(Default constructor)"));
2343 ExcMessage(
"Uninitialized TrilinosPayload::inv_Tvmult called "
2344 "(Default constructor)"));
2354 matrix.trilinos_matrix()),
2356 matrix_exemplar.trilinos_matrix().UseTranspose(),
2357 matrix_exemplar.get_mpi_communicator(),
2358 matrix_exemplar.locally_owned_domain_indices(),
2359 matrix_exemplar.locally_owned_range_indices())
2369 matrix.trilinos_matrix()),
2371 payload_exemplar.UseTranspose(),
2372 payload_exemplar.get_mpi_communicator(),
2373 payload_exemplar.locally_owned_domain_indices(),
2374 payload_exemplar.locally_owned_range_indices())
2384 matrix_exemplar.trilinos_matrix().UseTranspose(),
2385 matrix_exemplar.get_mpi_communicator(),
2386 matrix_exemplar.locally_owned_domain_indices(),
2387 matrix_exemplar.locally_owned_range_indices())
2396 preconditioner.trilinos_operator(),
2398 preconditioner_exemplar.trilinos_operator().UseTranspose(),
2399 preconditioner_exemplar.get_mpi_communicator(),
2400 preconditioner_exemplar.locally_owned_domain_indices(),
2401 preconditioner_exemplar.locally_owned_range_indices())
2411 payload_exemplar.UseTranspose(),
2412 payload_exemplar.get_mpi_communicator(),
2413 payload_exemplar.locally_owned_domain_indices(),
2414 payload_exemplar.locally_owned_range_indices())
2420 : vmult(payload.vmult)
2421 , Tvmult(payload.Tvmult)
2422 , inv_vmult(payload.inv_vmult)
2423 , inv_Tvmult(payload.inv_Tvmult)
2424 , use_transpose(payload.use_transpose)
2425 , communicator(payload.communicator)
2426 , domain_map(payload.domain_map)
2427 , range_map(payload.range_map)
2436 : use_transpose(false)
2439 communicator(first_op.communicator)
2440 , domain_map(second_op.domain_map)
2441 , range_map(first_op.range_map)
2452 tril_dst = tril_src;
2456 tril_dst = tril_src;
2460 tril_dst = tril_src;
2464 tril_dst = tril_src;
2478 const int ierr = tril_dst.PutScalar(0.0);
2484 const int ierr = tril_dst.PutScalar(0.0);
2491 ExcMessage(
"Cannot compute inverse of null operator"));
2493 const int ierr = tril_dst.PutScalar(0.0);
2500 ExcMessage(
"Cannot compute inverse of null operator"));
2502 const int ierr = tril_dst.PutScalar(0.0);
2602 return "TrilinosPayload";
2660 "Operators are set to work on incompatible IndexSets."));
2664 "Operators are set to work on incompatible IndexSets."));
2669 return_op.
vmult = [first_op, second_op](Range &tril_dst,
2670 const Domain &tril_src) {
2678 i->reinit(
IndexSet(first_op_init_map),
2683 const size_type i_local_size = i->end() - i->begin();
2687 Intermediate tril_int(View,
2695 second_op.
Apply(tril_src, tril_int);
2696 first_op.
Apply(tril_src, tril_dst);
2697 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2701 return_op.
Tvmult = [first_op, second_op](Domain &tril_dst,
2702 const Range &tril_src) {
2717 i->reinit(
IndexSet(first_op_init_map),
2722 const size_type i_local_size = i->end() - i->begin();
2726 Intermediate tril_int(View,
2734 second_op.
Apply(tril_src, tril_int);
2735 first_op.
Apply(tril_src, tril_dst);
2736 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2744 return_op.
inv_vmult = [first_op, second_op](Domain &tril_dst,
2745 const Range &tril_src) {
2753 i->reinit(
IndexSet(first_op_init_map),
2758 const size_type i_local_size = i->end() - i->begin();
2762 Intermediate tril_int(View,
2772 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2776 return_op.
inv_Tvmult = [first_op, second_op](Range &tril_dst,
2777 const Domain &tril_src) {
2792 i->reinit(
IndexSet(first_op_init_map),
2797 const size_type i_local_size = i->end() - i->begin();
2801 Intermediate tril_int(View,
2811 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2836 "Operators are set to work on incompatible IndexSets."));
2841 return_op.
vmult = [first_op, second_op](Range &tril_dst,
2842 const Domain &tril_src) {
2850 i->reinit(
IndexSet(first_op_init_map),
2855 const size_type i_local_size = i->end() - i->begin();
2859 Intermediate tril_int(View,
2867 second_op.
Apply(tril_src, tril_int);
2868 first_op.
Apply(tril_int, tril_dst);
2871 return_op.
Tvmult = [first_op, second_op](Domain &tril_dst,
2872 const Range &tril_src) {
2887 i->reinit(
IndexSet(first_op_init_map),
2892 const size_type i_local_size = i->end() - i->begin();
2896 Intermediate tril_int(View,
2903 first_op.
Apply(tril_src, tril_int);
2904 second_op.
Apply(tril_int, tril_dst);
2911 return_op.
inv_vmult = [first_op, second_op](Domain &tril_dst,
2912 const Range &tril_src) {
2920 i->reinit(
IndexSet(first_op_init_map),
2925 const size_type i_local_size = i->end() - i->begin();
2929 Intermediate tril_int(View,
2941 return_op.
inv_Tvmult = [first_op, second_op](Range &tril_dst,
2942 const Domain &tril_src) {
2957 i->reinit(
IndexSet(first_op_init_map),
2962 const size_type i_local_size = i->end() - i->begin();
2966 Intermediate tril_int(View,
2994# include "lac/trilinos_sparse_matrix.inst"
3009 const ::SparsityPattern &,
3025 const ::Vector<double> &)
const;
3030 const ::LinearAlgebra::distributed::Vector<double, MemorySpace::Host>
3036 const ::LinearAlgebra::distributed::Vector<float, MemorySpace::Host>
3039# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3040# if defined(HAVE_TPETRA_INST_DOUBLE)
3044 const ::LinearAlgebra::TpetraWrappers::Vector<
double,
3051 const ::LinearAlgebra::TpetraWrappers::Vector<
double,
3056# if defined(HAVE_TPETRA_INST_FLOAT)
3060 const ::LinearAlgebra::TpetraWrappers::Vector<
float,
3067 const ::LinearAlgebra::TpetraWrappers::Vector<
float,
3076 const ::LinearAlgebra::EpetraWrappers::Vector &)
const;
3083 const ::Vector<double> &)
const;
3088 const ::LinearAlgebra::distributed::Vector<double, MemorySpace::Host>
3091# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3092# if defined(HAVE_TPETRA_INST_DOUBLE)
3096 const ::LinearAlgebra::TpetraWrappers::Vector<
double,
3103 const ::LinearAlgebra::TpetraWrappers::Vector<
double,
3108# if defined(HAVE_TPETRA_INST_FLOAT)
3112 const ::LinearAlgebra::TpetraWrappers::Vector<
float,
3119 const ::LinearAlgebra::TpetraWrappers::Vector<
float,
3128 const ::LinearAlgebra::EpetraWrappers::Vector &)
const;
3135 const ::Vector<double> &)
const;
3140 const ::LinearAlgebra::distributed::Vector<double, MemorySpace::Host>
3143# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3144# if defined(HAVE_TPETRA_INST_DOUBLE)
3148 const ::LinearAlgebra::TpetraWrappers::Vector<
double,
3155 const ::LinearAlgebra::TpetraWrappers::Vector<
double,
3160# if defined(HAVE_TPETRA_INST_FLOAT)
3164 const ::LinearAlgebra::TpetraWrappers::Vector<
float,
3171 const ::LinearAlgebra::TpetraWrappers::Vector<
float,
3180 const ::LinearAlgebra::EpetraWrappers::Vector &)
const;
3187 const ::Vector<double> &)
const;
3192 const ::LinearAlgebra::distributed::Vector<double, MemorySpace::Host>
3198 const ::LinearAlgebra::distributed::Vector<float, MemorySpace::Host>
3202# ifdef DEAL_II_TRILINOS_WITH_TPETRA
3203# if defined(HAVE_TPETRA_INST_DOUBLE)
3207 const ::LinearAlgebra::TpetraWrappers::Vector<
double,
3214 const ::LinearAlgebra::TpetraWrappers::Vector<
double,
3219# if defined(HAVE_TPETRA_INST_FLOAT)
3223 const ::LinearAlgebra::TpetraWrappers::Vector<
float,
3230 const ::LinearAlgebra::TpetraWrappers::Vector<
float,
3239 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
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
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
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
void reinit(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const ::SparseMatrix< number > &dealii_sparse_matrix, const MPI_Comm communicator=MPI_COMM_WORLD, const double drop_tolerance=1e-13, const bool copy_values=true, const ::SparsityPattern *use_this_sparsity=nullptr)
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
::types::global_dof_index size_type
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_DISABLE_EXTRA_DIAGNOSTICS
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_NOT_IMPLEMENTED()
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)
IndexSet complete_index_set(const IndexSet::size_type N)
std::vector< index_type > data
@ 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)