19 #ifdef DEAL_II_WITH_TRILINOS 32 # include <boost/container/small_vector.hpp> 35 # ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT 36 # include <EpetraExt_MatrixMatrix.h> 38 # include <Epetra_Export.h> 39 # include <Teuchos_RCP.hpp> 40 # include <ml_epetra_utils.h> 41 # include <ml_struct.h> 51 template <
typename VectorType>
52 typename VectorType::value_type *
58 template <
typename VectorType>
59 const typename VectorType::value_type *
65 template <
typename VectorType>
66 typename VectorType::value_type *
72 template <
typename VectorType>
73 const typename VectorType::value_type *
79 # ifdef DEAL_II_WITH_MPI 108 # ifdef DEAL_II_TRILINOS_WITH_TPETRA 109 template <
typename Number>
116 template <
typename Number>
123 template <
typename Number>
131 template <
typename Number>
146 AccessorBase::visit_present_row()
154 if ((this->a_row ==
matrix->m()) ||
155 (
matrix->in_local_range(this->a_row) ==
false))
157 colnum_cache.reset();
166 if (value_cache.get() ==
nullptr)
169 std::make_shared<std::vector<TrilinosScalar>>(
matrix->n());
170 colnum_cache = std::make_shared<std::vector<size_type>>(
matrix->n());
174 value_cache->resize(
matrix->n());
175 colnum_cache->resize(
matrix->n());
178 int ierr =
matrix->trilinos_matrix().ExtractGlobalRowCopy(
184 colnum_cache->data()));
185 value_cache->resize(ncols);
186 colnum_cache->resize(ncols);
213 new Epetra_FECrsMatrix(
View, *column_space_map, *column_space_map, 0))
224 const unsigned int n_max_entries_per_row)
239 matrix(new Epetra_FECrsMatrix(
245 n_max_entries_per_row,
255 const std::vector<unsigned int> &n_entries_per_row)
260 ,
matrix(new Epetra_FECrsMatrix(
266 reinterpret_cast<
int *>(
267 const_cast<unsigned
int *>(n_entries_per_row.data())),
277 const unsigned int n_max_entries_per_row)
279 parallel_partitioning.make_trilinos_map(communicator, false)))
280 ,
matrix(new Epetra_FECrsMatrix(Copy,
282 n_max_entries_per_row,
292 const std::vector<unsigned int> &n_entries_per_row)
294 parallel_partitioning.make_trilinos_map(communicator, false)))
295 ,
matrix(new Epetra_FECrsMatrix(Copy,
297 reinterpret_cast<
int *>(
298 const_cast<unsigned
int *>(
299 n_entries_per_row.data())),
308 const IndexSet &col_parallel_partitioning,
312 col_parallel_partitioning.make_trilinos_map(communicator, false)))
313 ,
matrix(new Epetra_FECrsMatrix(
315 row_parallel_partitioning.make_trilinos_map(communicator, false),
316 n_max_entries_per_row,
325 const IndexSet &col_parallel_partitioning,
327 const std::vector<unsigned int> &n_entries_per_row)
329 col_parallel_partitioning.make_trilinos_map(communicator, false)))
330 ,
matrix(new Epetra_FECrsMatrix(
332 row_parallel_partitioning.make_trilinos_map(communicator, false),
333 reinterpret_cast<
int *>(
334 const_cast<unsigned
int *>(n_entries_per_row.data())),
345 new Epetra_FECrsMatrix(Copy,
353 "The Trilinos sparsity pattern has not been compressed."));
361 ,
matrix(std::move(other.matrix))
367 other.last_action = Zero;
368 other.compressed =
false;
385 bool needs_deep_copy =
390 if (!needs_deep_copy)
397 const int row_local =
matrix->RowMap().LID(
398 static_cast<TrilinosWrappers::types::int_type>(row));
401 int n_entries, rhs_n_entries;
403 int * index_ptr, *rhs_index_ptr;
404 int ierr = rhs.
matrix->ExtractMyRowView(row_local,
411 ierr =
matrix->ExtractMyRowView(row_local,
417 if (n_entries != rhs_n_entries ||
418 std::memcmp(static_cast<void *>(index_ptr),
419 static_cast<void *>(rhs_index_ptr),
420 sizeof(
int) * n_entries) != 0)
422 needs_deep_copy =
true;
426 for (
int i = 0; i < n_entries; ++i)
427 value_ptr[i] = rhs_value_ptr[i];
437 matrix = std::make_unique<Epetra_FECrsMatrix>(*rhs.
matrix);
444 std::make_unique<Epetra_CrsMatrix>(Copy, rhs.
nonlocal_matrix->Graph());
453 template <
typename SparsityPatternType>
455 reinit_matrix(
const IndexSet & row_parallel_partitioning,
456 const IndexSet & column_parallel_partitioning,
457 const SparsityPatternType & sparsity_pattern,
458 const bool exchange_data,
461 std::unique_ptr<Epetra_FECrsMatrix> &
matrix,
467 nonlocal_matrix.reset();
468 nonlocal_matrix_exporter.reset();
470 column_space_map = std::make_unique<Epetra_Map>(
473 if (column_space_map->Comm().MyPID() == 0)
476 row_parallel_partitioning.
size());
478 column_parallel_partitioning.
size());
481 Epetra_Map row_space_map =
491 trilinos_sparsity.
reinit(row_parallel_partitioning,
492 column_parallel_partitioning,
496 matrix = std::make_unique<Epetra_FECrsMatrix>(
505 std::vector<int> n_entries_per_row(last_row - first_row);
507 for (
size_type row = first_row; row < last_row; ++row)
508 n_entries_per_row[row - first_row] = sparsity_pattern.row_length(row);
524 std::unique_ptr<Epetra_CrsGraph> graph;
525 if (row_space_map.Comm().NumProc() > 1)
526 graph = std::make_unique<Epetra_CrsGraph>(Copy,
528 n_entries_per_row.data(),
531 graph = std::make_unique<Epetra_CrsGraph>(Copy,
534 n_entries_per_row.data(),
542 std::vector<TrilinosWrappers::types::int_type> row_indices;
544 for (
size_type row = first_row; row < last_row; ++row)
546 const int row_length = sparsity_pattern.row_length(row);
550 row_indices.resize(row_length, -1);
552 typename SparsityPatternType::iterator p =
553 sparsity_pattern.begin(row);
554 for (
size_type col = 0; p != sparsity_pattern.end(row); ++p, ++col)
555 row_indices[col] = p->column();
557 graph->Epetra_CrsGraph::InsertGlobalIndices(row,
566 graph->FillComplete(*column_space_map, row_space_map);
567 graph->OptimizeStorage();
575 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph,
false);
587 class Epetra_CrsGraphMod :
public Epetra_CrsGraph
590 Epetra_CrsGraphMod(
const Epetra_Map &row_map,
591 const int * n_entries_per_row)
592 : Epetra_CrsGraph(Copy, row_map, n_entries_per_row,
true)
596 SetIndicesAreGlobal()
598 this->Epetra_CrsGraph::SetIndicesAreGlobal(
true);
608 reinit_matrix(
const IndexSet & row_parallel_partitioning,
609 const IndexSet & column_parallel_partitioning,
611 const bool exchange_data,
613 std::unique_ptr<Epetra_Map> & column_space_map,
614 std::unique_ptr<Epetra_FECrsMatrix> &matrix,
615 std::unique_ptr<Epetra_CrsMatrix> & nonlocal_matrix,
616 std::unique_ptr<Epetra_Export> & nonlocal_matrix_exporter)
619 nonlocal_matrix.reset();
620 nonlocal_matrix_exporter.reset();
622 column_space_map = std::make_unique<Epetra_Map>(
626 row_parallel_partitioning.
size());
628 column_parallel_partitioning.
size());
630 Epetra_Map row_space_map =
635 if (relevant_rows.size() == 0)
639 relevant_rows.add_range(
642 relevant_rows.compress();
643 Assert(relevant_rows.n_elements() >=
644 static_cast<unsigned int>(row_space_map.NumMyElements()),
646 "Locally relevant rows of sparsity pattern must contain " 647 "all locally owned rows"));
652 const bool have_ghost_rows = [&]() {
653 std::vector<::types::global_dof_index> indices;
654 relevant_rows.fill_index_vector(indices);
655 Epetra_Map relevant_map(
663 row_space_map.Comm());
664 return !relevant_map.SameAs(row_space_map);
667 const unsigned int n_rows = relevant_rows.n_elements();
668 std::vector<TrilinosWrappers::types::int_type> ghost_rows;
669 std::vector<int> n_entries_per_row(row_space_map.NumMyElements());
670 std::vector<int> n_entries_per_ghost_row;
671 for (
unsigned int i = 0, own = 0; i < n_rows; ++i)
674 relevant_rows.nth_index_in_set(i);
675 if (row_space_map.MyGID(global_row))
676 n_entries_per_row[own++] = sparsity_pattern.
row_length(global_row);
677 else if (sparsity_pattern.
row_length(global_row) > 0)
679 ghost_rows.push_back(global_row);
680 n_entries_per_ghost_row.push_back(
685 Epetra_Map off_processor_map(-1,
687 (ghost_rows.size() > 0) ?
688 (ghost_rows.data()) :
691 row_space_map.Comm());
693 std::unique_ptr<Epetra_CrsGraph> graph;
694 std::unique_ptr<Epetra_CrsGraphMod> nonlocal_graph;
695 if (row_space_map.Comm().NumProc() > 1)
698 std::make_unique<Epetra_CrsGraph>(Copy,
700 (n_entries_per_row.size() > 0) ?
701 (n_entries_per_row.data()) :
703 exchange_data ?
false :
true);
704 if (have_ghost_rows ==
true)
705 nonlocal_graph = std::make_unique<Epetra_CrsGraphMod>(
706 off_processor_map, n_entries_per_ghost_row.data());
710 std::make_unique<Epetra_CrsGraph>(Copy,
713 (n_entries_per_row.size() > 0) ?
714 (n_entries_per_row.data()) :
719 std::vector<TrilinosWrappers::types::int_type> row_indices;
721 for (
unsigned int i = 0; i < n_rows; ++i)
724 relevant_rows.nth_index_in_set(i);
729 row_indices.resize(row_length, -1);
731 row_indices[col] = sparsity_pattern.
column_number(global_row, col);
733 if (row_space_map.MyGID(global_row))
734 graph->InsertGlobalIndices(global_row,
740 nonlocal_graph->InsertGlobalIndices(global_row,
747 if (nonlocal_graph.get() !=
nullptr)
753 nonlocal_graph->SetIndicesAreGlobal();
754 Assert(nonlocal_graph->IndicesAreGlobal() ==
true,
756 nonlocal_graph->FillComplete(*column_space_map, row_space_map);
757 nonlocal_graph->OptimizeStorage();
762 Epetra_Export exporter(nonlocal_graph->RowMap(), row_space_map);
763 int ierr = graph->Export(*nonlocal_graph, exporter, Add);
769 std::make_unique<Epetra_CrsMatrix>(Copy, *nonlocal_graph);
772 graph->FillComplete(*column_space_map, row_space_map);
773 graph->OptimizeStorage();
778 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph,
false);
784 template <
typename SparsityPatternType>
801 template <
typename SparsityPatternType>
802 inline typename std::enable_if<
803 !std::is_same<SparsityPatternType,
806 const IndexSet & col_parallel_partitioning,
807 const SparsityPatternType &sparsity_pattern,
809 const bool exchange_data)
811 reinit_matrix(row_parallel_partitioning,
812 col_parallel_partitioning,
838 matrix = std::make_unique<Epetra_FECrsMatrix>(
843 std::make_unique<Epetra_CrsMatrix>(Copy,
857 if (
this == &sparse_matrix)
861 std::make_unique<Epetra_Map>(sparse_matrix.
trilinos_matrix().DomainMap());
864 matrix = std::make_unique<Epetra_FECrsMatrix>(
879 template <
typename number>
882 const IndexSet & row_parallel_partitioning,
883 const IndexSet & col_parallel_partitioning,
884 const ::SparseMatrix<number> &dealii_sparse_matrix,
886 const double drop_tolerance,
887 const bool copy_values,
888 const ::SparsityPattern * use_this_sparsity)
890 if (copy_values ==
false)
894 if (use_this_sparsity ==
nullptr)
895 reinit(row_parallel_partitioning,
896 col_parallel_partitioning,
897 dealii_sparse_matrix.get_sparsity_pattern(),
901 reinit(row_parallel_partitioning,
902 col_parallel_partitioning,
909 const size_type n_rows = dealii_sparse_matrix.m();
914 const ::SparsityPattern &sparsity_pattern =
915 (use_this_sparsity !=
nullptr) ?
917 dealii_sparse_matrix.get_sparsity_pattern();
919 if (
matrix.get() ==
nullptr ||
m() != n_rows ||
922 reinit(row_parallel_partitioning,
923 col_parallel_partitioning,
935 std::vector<size_type> row_indices(maximum_row_length);
936 std::vector<TrilinosScalar>
values(maximum_row_length);
938 for (
size_type row = 0; row < n_rows; ++row)
940 if (row_parallel_partitioning.
is_element(row) ==
true)
943 sparsity_pattern.begin(row);
944 typename ::SparseMatrix<number>::const_iterator it =
945 dealii_sparse_matrix.begin(row);
947 if (sparsity_pattern.n_rows() == sparsity_pattern.n_cols())
951 if (
std::fabs(it->value()) > drop_tolerance)
953 values[col] = it->value();
954 row_indices[col++] = it->column();
960 while (it != dealii_sparse_matrix.end(row) &&
961 select_index != sparsity_pattern.end(row))
963 while (select_index->column() < it->column() &&
964 select_index != sparsity_pattern.end(row))
966 while (it->column() < select_index->column() &&
967 it != dealii_sparse_matrix.end(row))
970 if (it == dealii_sparse_matrix.end(row))
972 if (
std::fabs(it->value()) > drop_tolerance)
974 values[col] = it->value();
975 row_indices[col++] = it->column();
982 reinterpret_cast<size_type *
>(row_indices.data()),
991 template <
typename number>
994 const ::SparseMatrix<number> &dealii_sparse_matrix,
995 const double drop_tolerance,
996 const bool copy_values,
997 const ::SparsityPattern * use_this_sparsity)
1001 dealii_sparse_matrix,
1012 const bool copy_values)
1014 Assert(input_matrix.Filled() ==
true,
1015 ExcMessage(
"Input CrsMatrix has not called FillComplete()!"));
1019 const Epetra_CrsGraph *graph = &input_matrix.Graph();
1024 matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph,
false);
1026 matrix->FillComplete(*column_space_map, input_matrix.RangeMap(),
true);
1028 if (copy_values ==
true)
1034 const size_type my_nonzeros = input_matrix.NumMyNonzeros();
1035 std::memcpy(values, in_values, my_nonzeros *
sizeof(
TrilinosScalar));
1059 "compress() can only be called with VectorOperation add, insert, or unknown"));
1067 ExcMessage(
"Operation and argument to compress() do not match"));
1093 ierr =
matrix->OptimizeStorage();
1115 matrix->FillComplete();
1131 matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
1138 matrix->ExtractMyRowView(local_row, num_entries, values, col_indices);
1143 const std::ptrdiff_t diag_index =
1144 std::find(col_indices, col_indices + num_entries, local_row) -
1148 if (diag_index != j || new_diag_value == 0)
1151 if (diag_index != num_entries)
1152 values[diag_index] = new_diag_value;
1162 for (
const auto row : rows)
1174 matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
1176 matrix->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
1185 if (trilinos_i == -1)
1200 int nnz_present =
matrix->NumMyEntries(trilinos_i);
1209 int ierr =
matrix->ExtractMyRowView(trilinos_i,
1216 Assert(nnz_present == nnz_extracted,
1222 const std::ptrdiff_t local_col_index =
1223 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1234 if (local_col_index == nnz_present)
1239 value = values[local_col_index];
1253 matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
1255 matrix->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
1265 if ((trilinos_i == -1) || (trilinos_j == -1))
1276 int nnz_present =
matrix->NumMyEntries(trilinos_i);
1284 int ierr =
matrix->ExtractMyRowView(trilinos_i,
1291 Assert(nnz_present == nnz_extracted,
1297 const std::ptrdiff_t local_col_index =
1298 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1308 if (local_col_index == nnz_present)
1311 value = values[local_col_index];
1349 matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
1357 int ierr =
matrix->NumMyRowEntries(local_row, ncols);
1361 return static_cast<unsigned int>(ncols);
1368 const std::vector<size_type> & col_indices,
1370 const bool elide_zero_values)
1372 Assert(row_indices.size() == values.
m(),
1374 Assert(col_indices.size() == values.
n(),
1377 for (
size_type i = 0; i < row_indices.size(); ++i)
1389 const std::vector<size_type> & col_indices,
1390 const std::vector<TrilinosScalar> &
values,
1391 const bool elide_zero_values)
1393 Assert(col_indices.size() == values.size(),
1411 const bool elide_zero_values)
1431 boost::container::small_vector<TrilinosScalar, 200> local_value_array(
1433 boost::container::small_vector<TrilinosWrappers::types::int_type, 200>
1434 local_index_array(n_cols);
1439 if (elide_zero_values ==
false)
1451 col_index_ptr = local_index_array.data();
1452 col_value_ptr = local_value_array.data();
1457 const double value = values[j];
1461 local_index_array[n_columns] = col_indices[j];
1462 local_value_array[n_columns] = value;
1479 if (
matrix->RowMap().MyGID(
1480 static_cast<TrilinosWrappers::types::int_type>(row)) ==
true)
1482 if (
matrix->Filled() ==
false)
1484 ierr =
matrix->Epetra_CrsMatrix::InsertGlobalValues(
1485 row, static_cast<int>(n_columns), col_value_ptr, col_index_ptr);
1494 ierr =
matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row,
1509 if (
matrix->Filled() ==
false)
1511 ierr =
matrix->InsertGlobalValues(1,
1516 Epetra_FECrsMatrix::ROW_MAJOR);
1521 ierr =
matrix->ReplaceGlobalValues(1,
1526 Epetra_FECrsMatrix::ROW_MAJOR);
1543 const bool elide_zero_values)
1545 Assert(indices.size() == values.
m(),
1549 for (
size_type i = 0; i < indices.size(); ++i)
1561 const std::vector<size_type> & col_indices,
1563 const bool elide_zero_values)
1565 Assert(row_indices.size() == values.
m(),
1567 Assert(col_indices.size() == values.
n(),
1570 for (
size_type i = 0; i < row_indices.size(); ++i)
1582 const std::vector<size_type> & col_indices,
1583 const std::vector<TrilinosScalar> &values,
1584 const bool elide_zero_values)
1586 Assert(col_indices.size() == values.size(),
1603 const bool elide_zero_values,
1625 boost::container::small_vector<TrilinosScalar, 100> local_value_array(
1627 boost::container::small_vector<TrilinosWrappers::types::int_type, 100>
1628 local_index_array(n_cols);
1633 if (elide_zero_values ==
false)
1649 col_index_ptr = local_index_array.data();
1650 col_value_ptr = local_value_array.data();
1655 const double value = values[j];
1660 local_index_array[n_columns] = col_indices[j];
1661 local_value_array[n_columns] = value;
1677 if (
matrix->RowMap().MyGID(
1678 static_cast<TrilinosWrappers::types::int_type>(row)) ==
true)
1680 ierr =
matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row,
1693 static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1694 ExcMessage(
"Attempted to write into off-processor matrix row " 1695 "that has not be specified as being writable upon " 1713 ierr =
matrix->SumIntoGlobalValues(1,
1718 Epetra_FECrsMatrix::ROW_MAJOR);
1725 std::cout <<
"------------------------------------------" << std::endl;
1726 std::cout <<
"Got error " << ierr <<
" in row " << row <<
" of proc " 1727 <<
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(
1734 static_cast<TrilinosWrappers::types::int_type>(row)) ==
1738 <<
" has the following indices:" << std::endl;
1739 std::vector<TrilinosWrappers::types::int_type> indices;
1740 const Epetra_CrsGraph * graph =
1743 static_cast<TrilinosWrappers::types::int_type>(row)) ==
false) ?
1747 indices.resize(graph->NumGlobalIndices(row));
1749 graph->ExtractGlobalRowCopy(row,
1756 std::cout << indices[i] <<
" ";
1757 std::cout << std::endl << std::endl;
1774 const int ierr =
matrix->PutScalar(d);
1792 ExcMessage(
"Can only add matrices with same distribution of rows"));
1794 ExcMessage(
"Addition of matrices only allowed if matrices are " 1795 "filled, i.e., compress() has been called"));
1797 const bool same_col_map =
matrix->ColMap().SameAs(rhs.
matrix->ColMap());
1801 const int row_local =
matrix->RowMap().LID(
1802 static_cast<TrilinosWrappers::types::int_type>(row));
1809 int n_entries, rhs_n_entries;
1811 int * index_ptr, *rhs_index_ptr;
1812 int ierr = rhs.
matrix->ExtractMyRowView(row_local,
1820 matrix->ExtractMyRowView(row_local, n_entries, value_ptr, index_ptr);
1822 bool expensive_checks = (n_entries != rhs_n_entries || !same_col_map);
1823 if (!expensive_checks)
1827 expensive_checks = std::memcmp(static_cast<void *>(index_ptr),
1828 static_cast<void *>(rhs_index_ptr),
1829 sizeof(
int) * n_entries) != 0;
1830 if (!expensive_checks)
1831 for (
int i = 0; i < n_entries; ++i)
1832 value_ptr[i] += rhs_value_ptr[i] * factor;
1838 if (expensive_checks)
1840 for (
int i = 0; i < rhs_n_entries; ++i)
1842 if (rhs_value_ptr[i] == 0.)
1846 int local_col =
matrix->ColMap().LID(rhs_global_col);
1848 index_ptr + n_entries,
1850 Assert(local_index != index_ptr + n_entries &&
1851 *local_index == local_col,
1853 "Adding the entries from the other matrix " 1854 "failed, because the sparsity pattern " 1855 "of that matrix includes more elements than the " 1856 "calling matrix, which is not allowed."));
1857 value_ptr[local_index - index_ptr] += factor * rhs_value_ptr[i];
1875 if (!
matrix->UseTranspose())
1877 ierr =
matrix->SetUseTranspose(
true);
1882 ierr =
matrix->SetUseTranspose(
false);
1892 const int ierr =
matrix->Scale(a);
1908 const int ierr =
matrix->Scale(factor);
1921 return matrix->NormOne();
1930 return matrix->NormInf();
1939 return matrix->NormFrobenius();
1946 namespace SparseMatrixImplementation
1948 template <
typename VectorType>
1962 "Column map of matrix does not fit with vector map!"));
1964 ExcMessage(
"Row map of matrix does not fit with vector map!"));
1973 template <
typename VectorType>
1974 typename std::enable_if<
1975 std::is_same<typename VectorType::value_type, TrilinosScalar>::value>::type
1991 Epetra_MultiVector tril_dst(
1993 Epetra_MultiVector tril_src(
View,
2000 const int ierr =
matrix->Multiply(
false, tril_src, tril_dst);
2007 template <
typename VectorType>
2008 typename std::enable_if<
2009 !std::is_same<typename VectorType::value_type, TrilinosScalar>::value>::type
2017 template <
typename VectorType>
2018 typename std::enable_if<
2019 std::is_same<typename VectorType::value_type, TrilinosScalar>::value>::type
2033 Epetra_MultiVector tril_dst(
2035 Epetra_MultiVector tril_src(
View,
2041 const int ierr =
matrix->Multiply(
true, tril_src, tril_dst);
2048 template <
typename VectorType>
2049 typename std::enable_if<
2050 !std::is_same<typename VectorType::value_type, TrilinosScalar>::value>::type
2058 template <
typename VectorType>
2067 tmp_vector.reinit(dst,
true);
2068 vmult(tmp_vector, src);
2074 template <
typename VectorType>
2083 tmp_vector.reinit(dst,
true);
2096 temp_vector.
reinit(v,
true);
2098 vmult(temp_vector, v);
2099 return temp_vector * v;
2111 temp_vector.
reinit(v,
true);
2113 vmult(temp_vector, v);
2114 return u * temp_vector;
2142 const bool transpose_left)
2144 const bool use_vector = (V.
size() == inputright.
m() ? true :
false);
2145 if (transpose_left ==
false)
2147 Assert(inputleft.
n() == inputright.
m(),
2151 ExcMessage(
"Parallel partitioning of A and B does not fit."));
2155 Assert(inputleft.
m() == inputright.
m(),
2159 ExcMessage(
"Parallel partitioning of A and B does not fit."));
2170 Teuchos::RCP<Epetra_CrsMatrix> mod_B;
2171 if (use_vector ==
false)
2173 mod_B = Teuchos::rcp(const_cast<Epetra_CrsMatrix *>(
2179 mod_B = Teuchos::rcp(
2185 ExcMessage(
"Parallel distribution of matrix B and vector V " 2186 "does not match."));
2189 for (
int i = 0; i < local_N; ++i)
2192 double *new_data, *B_data;
2193 mod_B->ExtractMyRowView(i, N_entries, new_data);
2199 new_data[j] = value * B_data[j];
2210 # ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT 2215 const_cast<Epetra_CrsMatrix &
>(
2216 tmp_result.trilinos_matrix()));
2220 result.
reinit(tmp_result.trilinos_matrix());
2258 const bool print_detailed_trilinos_information)
const 2260 if (print_detailed_trilinos_information ==
true)
2268 for (
int i = 0; i < matrix->NumMyRows(); ++i)
2271 matrix->ExtractMyRowView(i, num_entries, values, indices);
2279 <<
") " << values[j] << std::endl;
2292 sizeof(*this) +
sizeof(*matrix) +
sizeof(*
matrix->Graph().DataPtr());
2295 matrix->NumMyNonzeros() +
2304 # ifdef DEAL_II_WITH_MPI 2306 const Epetra_MpiComm *mpi_comm =
2307 dynamic_cast<const Epetra_MpiComm *
>(&
matrix->RangeMap().Comm());
2309 return mpi_comm->Comm();
2312 return MPI_COMM_SELF;
2325 # ifndef DEAL_II_WITH_MPI 2327 make_serial_Epetra_map(
const IndexSet &serial_partitioning)
2334 Epetra_SerialComm());
2339 namespace LinearOperatorImplementation
2341 TrilinosPayload::TrilinosPayload()
2342 : use_transpose(false)
2344 # ifdef DEAL_II_WITH_MPI
2345 communicator(MPI_COMM_SELF)
2346 , domain_map(
IndexSet().make_trilinos_map(communicator.Comm()))
2347 , range_map(
IndexSet().make_trilinos_map(communicator.Comm()))
2355 ExcMessage(
"Uninitialized TrilinosPayload::vmult called " 2356 "(Default constructor)"));
2361 ExcMessage(
"Uninitialized TrilinosPayload::Tvmult called " 2362 "(Default constructor)"));
2367 ExcMessage(
"Uninitialized TrilinosPayload::inv_vmult called " 2368 "(Default constructor)"));
2373 ExcMessage(
"Uninitialized TrilinosPayload::inv_Tvmult called " 2374 "(Default constructor)"));
2385 # ifdef DEAL_II_WITH_MPI
2401 const Domain &tril_src) {
2403 Assert(&tril_src != &tril_dst,
2423 const Range &tril_src) {
2425 Assert(&tril_src != &tril_dst,
2440 Epetra_CrsMatrix &tril_mtrx_non_const =
2442 tril_mtrx_non_const.SetUseTranspose(
2446 tril_mtrx_non_const.SetUseTranspose(
2452 ExcMessage(
"Uninitialized TrilinosPayload::inv_vmult called " 2453 "(Matrix constructor with matrix exemplar)"));
2458 ExcMessage(
"Uninitialized TrilinosPayload::inv_Tvmult called " 2459 "(Matrix constructor with matrix exemplar)"));
2470 # ifdef DEAL_II_WITH_MPI
2485 vmult = [&matrix_exemplar, &preconditioner](
Range & tril_dst,
2486 const Domain &tril_src) {
2489 Assert(&tril_src != &tril_dst,
2507 Tvmult = [&matrix_exemplar, &preconditioner](
Domain & tril_dst,
2508 const Range &tril_src) {
2511 Assert(&tril_src != &tril_dst,
2534 const Range &tril_src) {
2537 Assert(&tril_src != &tril_dst,
2556 &preconditioner](
Range & tril_dst,
2557 const Domain &tril_src) {
2560 Assert(&tril_src != &tril_dst,
2589 preconditioner_exemplar.trilinos_operator().
UseTranspose())
2591 # ifdef DEAL_II_WITH_MPI
2604 vmult = [&preconditioner_exemplar,
2605 &preconditioner](
Range &tril_dst,
const Domain &tril_src) {
2608 Assert(&tril_src != &tril_dst,
2626 Tvmult = [&preconditioner_exemplar,
2627 &preconditioner](
Domain &tril_dst,
const Range &tril_src) {
2630 Assert(&tril_src != &tril_dst,
2653 &preconditioner](
Domain &tril_dst,
const Range &tril_src) {
2656 Assert(&tril_src != &tril_dst,
2675 &preconditioner](
Range & tril_dst,
2676 const Domain &tril_src) {
2679 Assert(&tril_src != &tril_dst,
2737 tril_dst = tril_src;
2741 tril_dst = tril_src;
2745 tril_dst = tril_src;
2749 tril_dst = tril_src;
2763 const int ierr = tril_dst.PutScalar(0.0);
2769 const int ierr = tril_dst.PutScalar(0.0);
2776 ExcMessage(
"Cannot compute inverse of null operator"));
2778 const int ierr = tril_dst.PutScalar(0.0);
2785 ExcMessage(
"Cannot compute inverse of null operator"));
2787 const int ierr = tril_dst.PutScalar(0.0);
2826 # ifdef DEAL_II_WITH_MPI 2829 return MPI_COMM_SELF;
2891 return "TrilinosPayload";
2949 "Operators are set to work on incompatible IndexSets."));
2953 "Operators are set to work on incompatible IndexSets."));
2958 return_op.vmult = [first_op, second_op](
Range & tril_dst,
2959 const Domain &tril_src) {
2967 i->reinit(
IndexSet(first_op_init_map),
2972 const size_type i_local_size = i->end() - i->begin();
2976 (void)second_op_init_map;
2977 Intermediate tril_int(
View,
2979 const_cast<TrilinosScalar *>(i->begin()),
2985 second_op.
Apply(tril_src, tril_int);
2986 first_op.
Apply(tril_src, tril_dst);
2987 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2991 return_op.Tvmult = [first_op, second_op](
Domain & tril_dst,
2992 const Range &tril_src) {
3007 i->reinit(
IndexSet(first_op_init_map),
3012 const size_type i_local_size = i->end() - i->begin();
3016 (void)second_op_init_map;
3017 Intermediate tril_int(
View,
3019 const_cast<TrilinosScalar *>(i->begin()),
3025 second_op.
Apply(tril_src, tril_int);
3026 first_op.
Apply(tril_src, tril_dst);
3027 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
3035 return_op.inv_vmult = [first_op, second_op](
Domain & tril_dst,
3036 const Range &tril_src) {
3044 i->reinit(
IndexSet(first_op_init_map),
3049 const size_type i_local_size = i->end() - i->begin();
3053 (void)second_op_init_map;
3054 Intermediate tril_int(
View,
3056 const_cast<TrilinosScalar *>(i->begin()),
3064 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
3068 return_op.inv_Tvmult = [first_op, second_op](
Range & tril_dst,
3069 const Domain &tril_src) {
3084 i->reinit(
IndexSet(first_op_init_map),
3089 const size_type i_local_size = i->end() - i->begin();
3093 (void)second_op_init_map;
3094 Intermediate tril_int(
View,
3096 const_cast<TrilinosScalar *>(i->begin()),
3104 const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
3128 "Operators are set to work on incompatible IndexSets."));
3133 return_op.vmult = [first_op, second_op](
Range & tril_dst,
3134 const Domain &tril_src) {
3142 i->reinit(
IndexSet(first_op_init_map),
3147 const size_type i_local_size = i->end() - i->begin();
3151 (void)second_op_init_map;
3152 Intermediate tril_int(
View,
3154 const_cast<TrilinosScalar *>(i->begin()),
3160 second_op.
Apply(tril_src, tril_int);
3161 first_op.
Apply(tril_int, tril_dst);
3164 return_op.Tvmult = [first_op, second_op](
Domain & tril_dst,
3165 const Range &tril_src) {
3180 i->reinit(
IndexSet(first_op_init_map),
3185 const size_type i_local_size = i->end() - i->begin();
3189 (void)second_op_init_map;
3190 Intermediate tril_int(
View,
3192 const_cast<TrilinosScalar *>(i->begin()),
3197 first_op.
Apply(tril_src, tril_int);
3198 second_op.
Apply(tril_int, tril_dst);
3205 return_op.inv_vmult = [first_op, second_op](
Domain & tril_dst,
3206 const Range &tril_src) {
3214 i->reinit(
IndexSet(first_op_init_map),
3219 const size_type i_local_size = i->end() - i->begin();
3223 (void)second_op_init_map;
3224 Intermediate tril_int(
View,
3226 const_cast<TrilinosScalar *>(i->begin()),
3236 return_op.inv_Tvmult = [first_op, second_op](
Range & tril_dst,
3237 const Domain &tril_src) {
3252 i->reinit(
IndexSet(first_op_init_map),
3257 const size_type i_local_size = i->end() - i->begin();
3261 (void)second_op_init_map;
3262 Intermediate tril_int(
View,
3264 const_cast<TrilinosScalar *>(i->begin()),
3290 # include "trilinos_sparse_matrix.inst" 3305 const ::SparsityPattern &,
3321 const ::Vector<double> &)
const;
3326 const ::LinearAlgebra::distributed::Vector<double> &)
const;
3328 # ifdef DEAL_II_WITH_MPI 3329 # ifdef DEAL_II_TRILINOS_WITH_TPETRA 3333 const ::LinearAlgebra::TpetraWrappers::Vector<double> &)
const;
3338 const ::LinearAlgebra::TpetraWrappers::Vector<float> &)
const;
3352 const ::Vector<double> &)
const;
3357 const ::LinearAlgebra::distributed::Vector<double> &)
const;
3359 # ifdef DEAL_II_WITH_MPI 3360 # ifdef DEAL_II_TRILINOS_WITH_TPETRA 3364 const ::LinearAlgebra::TpetraWrappers::Vector<double> &)
const;
3369 const ::LinearAlgebra::TpetraWrappers::Vector<float> &)
const;
3383 const ::Vector<double> &)
const;
3388 const ::LinearAlgebra::distributed::Vector<double> &)
const;
3390 # ifdef DEAL_II_WITH_MPI 3391 # ifdef DEAL_II_TRILINOS_WITH_TPETRA 3395 const ::LinearAlgebra::TpetraWrappers::Vector<double> &)
const;
3400 const ::LinearAlgebra::TpetraWrappers::Vector<float> &)
const;
3414 const ::Vector<double> &)
const;
3419 const ::LinearAlgebra::distributed::Vector<double> &)
const;
3421 # ifdef DEAL_II_WITH_MPI 3422 # ifdef DEAL_II_TRILINOS_WITH_TPETRA 3426 const ::LinearAlgebra::TpetraWrappers::Vector<double> &)
const;
3431 const ::LinearAlgebra::TpetraWrappers::Vector<float> &)
const;
3444 #endif // DEAL_II_WITH_TRILINOS TrilinosScalar diag_element(const size_type i) const
Iterator lower_bound(Iterator first, Iterator last, const T &val)
static ::ExceptionBase & ExcTrilinosError(int arg1)
TrilinosScalar el(const size_type i, const size_type j) const
virtual double NormInf() const override
const Epetra_Map & domain_partitioner() const
TrilinosScalar frobenius_norm() const
void mmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
real_type l2_norm() const
#define AssertDimension(dim1, dim2)
TrilinosPayload null_payload() const
typename ::internal::FEValuesViews::ViewType< dim, spacedim, Extractor >::type View
virtual int SetUseTranspose(bool UseTranspose) override
Contents is actually a matrix.
static ::ExceptionBase & ExcSourceEqualsDestination()
TrilinosPayload operator*(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
void check_vector_map_equality(const Epetra_CrsMatrix &, const VectorType &, const VectorType &)
static ::ExceptionBase & ExcIO()
void clear_row(const size_type row, const TrilinosScalar new_diag_value=0)
size_type column_number(const size_type row, const size_type index) const
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
void vmult_add(VectorType &dst, const VectorType &src) const
Epetra_MultiVector VectorType
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
::types::global_dof_index size_type
Epetra_CombineMode last_action
const Epetra_FECrsGraph & trilinos_sparsity_pattern() const
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
#define AssertIndexRange(index, range)
TrilinosWrappers::types::int_type min_my_gid(const Epetra_BlockMap &map)
std::pair< size_type, size_type > local_range() const
virtual const char * Label() const override
TrilinosScalar operator()(const size_type i, const size_type j) const
std::function< void(VectorType &, const VectorType &)> Tvmult
TrilinosPayload operator+(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define AssertThrow(cond, exc)
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
TrilinosScalar residual(MPI::Vector &dst, const MPI::Vector &x, const MPI::Vector &b) const
const Epetra_Comm & comm_self()
void check_vector_map_equality(const Epetra_CrsMatrix &m, const TrilinosWrappers::MPI::Vector &in, const TrilinosWrappers::MPI::Vector &out)
static ::ExceptionBase & ExcTrilinosError(int arg1)
TrilinosScalar matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
static ::ExceptionBase & ExcDivideByZero()
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
void clear_rows(const std::vector< size_type > &rows, const TrilinosScalar new_diag_value=0)
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
SparseMatrix & operator*=(const TrilinosScalar factor)
std::unique_ptr< Epetra_FECrsMatrix > matrix
virtual int Apply(const VectorType &X, VectorType &Y) const override
SparseMatrix & operator=(const SparseMatrix &)=delete
TrilinosPayload transpose_payload() const
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1)
::types::global_dof_index size_type
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual const Epetra_Map & OperatorRangeMap() const override
TrilinosPayload identity_payload() const
void add(const size_type i, const size_type j, const TrilinosScalar value)
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::unique_ptr< Epetra_Export > nonlocal_matrix_exporter
TrilinosWrappers::types::int_type n_global_cols(const Epetra_CrsGraph &graph)
size_type memory_consumption() const
TrilinosWrappers::types::int_type global_row_index(const Epetra_CrsMatrix &matrix, const ::types::global_dof_index i)
const IndexSet & row_index_set() const
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
Epetra_MpiComm communicator
#define DEAL_II_NAMESPACE_CLOSE
std::function< void(VectorType &, const VectorType &)> inv_vmult
VectorType::value_type * end(VectorType &V)
TrilinosScalar matrix_norm_square(const MPI::Vector &v) const
const Epetra_FEVector & trilinos_vector() const
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
IndexSet locally_owned_range_indices() const
Expression fabs(const Expression &x)
const Epetra_CrsMatrix & trilinos_matrix() const
virtual bool UseTranspose() const override
TrilinosWrappers::types::int_type global_column_index(const Epetra_CrsMatrix &matrix, const ::types::global_dof_index i)
Epetra_Operator & trilinos_operator() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
IndexSet locally_owned_range_indices() const
virtual const Epetra_Map & OperatorDomainMap() const override
size_type n_nonzero_elements() const
std::enable_if< std::is_same< typename VectorType::value_type, TrilinosScalar >::value >::type Tvmult(VectorType &dst, const VectorType &src) const
MPI_Comm get_mpi_communicator() const
static ::ExceptionBase & ExcNotQuadratic()
TrilinosWrappers::types::int_type n_global_elements(const Epetra_BlockMap &map)
unsigned int row_length(const size_type row) const
void set_size(const size_type size)
IndexSet complete_index_set(const IndexSet::size_type N)
void Tvmult_add(VectorType &dst, const VectorType &src) const
unsigned int global_dof_index
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
std::enable_if< std::is_same< typename VectorType::value_type, TrilinosScalar >::value >::type vmult(VectorType &dst, const VectorType &src) const
const Tpetra::Vector< Number, int, types::global_dof_index > & trilinos_vector() const
void reinit(const SparsityPatternType &sparsity_pattern)
size_type row_length(const size_type row) const
TrilinosScalar linfty_norm() const
const Epetra_BlockMap & trilinos_partitioner() const
void copy_from(const SparseMatrix &source)
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_NAMESPACE_OPEN
VectorType::value_type * begin(VectorType &V)
std::function< void(VectorType &, const VectorType &)> vmult
SparseMatrix & operator/=(const TrilinosScalar factor)
TrilinosWrappers::types::int_type max_my_gid(const Epetra_BlockMap &map)
virtual const Epetra_Comm & Comm() const override
std::pair< size_type, size_type > local_range() const
static ::ExceptionBase & ExcNotImplemented()
IndexSet locally_owned_domain_indices() const
static ::ExceptionBase & ExcMatrixNotCompressed()
IndexSet locally_owned_domain_indices() const
const Epetra_CrsGraph & trilinos_sparsity_pattern() const
bool is_element(const size_type index) const
void Tmmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
std::unique_ptr< Epetra_Map > column_space_map
void perform_mmult(const SparseMatrix &inputleft, const SparseMatrix &inputright, SparseMatrix &result, const MPI::Vector &V, const bool transpose_left)
virtual bool HasNormInf() const override
void compress(::VectorOperation::values operation)
const Epetra_MultiVector & trilinos_vector() const
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
MPI_Comm get_mpi_communicator() const
void reinit(const size_type m, const size_type n, const size_type n_entries_per_row=0)
void set(const size_type i, const size_type j, const TrilinosScalar value)
size_type n_elements() const
#define AssertIsFinite(number)
unsigned int local_size() const
TrilinosScalar l1_norm() const
void check_vector_map_equality(const Epetra_CrsMatrix &mtrx, const Epetra_MultiVector &src, const Epetra_MultiVector &dst, const bool transpose)
static ::ExceptionBase & ExcInternalError()