From ec280f30600e81953ff12ee32aca46094a4ec77c Mon Sep 17 00:00:00 2001 From: Sebastian Kinnewig Date: Sun, 18 Feb 2024 15:10:04 +0100 Subject: [PATCH] Rework reinit_matrix in LA:TpetraWrappers::SparseMatrix --- .../trilinos_tpetra_sparse_matrix.templates.h | 130 ++++++++++++++++-- 1 file changed, 120 insertions(+), 10 deletions(-) diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h index 50deeda6de..650c1feb20 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h @@ -67,6 +67,118 @@ namespace LinearAlgebra // release memory before reallocation matrix.reset(); + // Get the Tpetra::Maps + Teuchos::RCP> row_space_map = + row_parallel_partitioning.make_tpetra_map_rcp(communicator, false); + + column_space_map = + column_parallel_partitioning.make_tpetra_map_rcp(communicator, false); + + if (column_space_map->getComm()->getRank() == 0) + { + AssertDimension(sparsity_pattern.n_rows(), + row_parallel_partitioning.size()); + AssertDimension(sparsity_pattern.n_cols(), + column_parallel_partitioning.size()); + } + + // if we want to exchange data, build a usual Trilinos sparsity pattern + // and let that handle the exchange. otherwise, manually create a + // CrsGraph, which consumes considerably less memory because it can set + // correct number of indices right from the start + if (exchange_data) + { + SparsityPattern trilinos_sparsity; + trilinos_sparsity.reinit(row_parallel_partitioning, + column_parallel_partitioning, + sparsity_pattern, + communicator, + exchange_data); + matrix = Utilities::Trilinos::internal::make_rcp< + MatrixType>( + trilinos_sparsity.trilinos_sparsity_pattern()); + + return; + } + + // compute the number of entries per row + const size_type first_row = row_space_map->getMinGlobalIndex(); + const size_type last_row = row_space_map->getMaxGlobalIndex() + 1; + + Teuchos::Array n_entries_per_row(last_row - first_row); + for (size_type row = first_row; row < last_row; ++row) + n_entries_per_row[row - first_row] = sparsity_pattern.row_length(row); + + // The deal.II notation of a Sparsity pattern corresponds to the + // Tpetra concept of a Graph. Hence, we generate a graph by copying + // the sparsity pattern into it, and then build up the matrix from the + // graph. This is considerable faster than directly filling elements + // into the matrix. Moreover, it consumes less memory, since the + // internal reordering is done on ints only, and we can leave the + // doubles aside. +# if DEAL_II_TRILINOS_VERSION_GTE(12, 16, 0) + Teuchos::RCP> graph = + Utilities::Trilinos::internal::make_rcp>( + row_space_map, n_entries_per_row); +# else + Teuchos::RCP> graph = + Utilities::Trilinos::internal::make_rcp>( + row_space_map, Teuchos::arcpFromArray(n_entries_per_row)); +# endif + + // This functions assumes that the sparsity pattern sits on all + // processors (completely). The parallel version uses a Tpetra graph + // that is already distributed. + + // now insert the indices + Teuchos::Array row_indices; + + for (size_type global_row = first_row; global_row < last_row; + ++global_row) + { + const int row_length = sparsity_pattern.row_length(global_row); + if (row_length == 0) + continue; + + row_indices.resize(row_length, -1); + for (size_type col = 0; col < row_length; ++col) + row_indices[col] = + sparsity_pattern.column_number(global_row, col); + + AssertIndexRange(global_row, row_space_map->getGlobalNumElements()); + graph->insertGlobalIndices(global_row, row_indices); + } + + // Eventually, optimize the graph structure (sort indices, make memory + // contiguous, etc.). note that the documentation of the function indeed + // states that we first need to provide the column (domain) map and then + // the row (range) map + graph->fillComplete(column_space_map, row_space_map); + + // check whether we got the number of columns right. + AssertDimension(sparsity_pattern.n_cols(), graph->getGlobalNumCols()); + + // And now finally generate the matrix. + matrix = + Utilities::Trilinos::internal::make_rcp>( + graph); + } + + + + template + void + reinit_matrix(const IndexSet &row_parallel_partitioning, + const IndexSet &column_parallel_partitioning, + const DynamicSparsityPattern &sparsity_pattern, + const bool exchange_data, + const MPI_Comm communicator, + Teuchos::RCP> &column_space_map, + Teuchos::RCP> &matrix) + { + // release memory before reallocation + matrix.reset(); + // Get the Tpetra::Maps Teuchos::RCP> row_space_map = row_parallel_partitioning.make_tpetra_map_rcp(communicator, false); @@ -135,14 +247,14 @@ namespace LinearAlgebra // into the matrix. Moreover, it consumes less memory, since the // internal reordering is done on ints only, and we can leave the // doubles aside. - Teuchos::RCP> graph; - # if DEAL_II_TRILINOS_VERSION_GTE(12, 16, 0) - graph = Utilities::Trilinos::internal::make_rcp>( - row_space_map, n_entries_per_row); + Teuchos::RCP> graph = + Utilities::Trilinos::internal::make_rcp>( + row_space_map, n_entries_per_row); # else - graph = Utilities::Trilinos::internal::make_rcp>( - row_space_map, Teuchos::arcpFromArray(n_entries_per_row)); + Teuchos::RCP> graph = + Utilities::Trilinos::internal::make_rcp>( + row_space_map, Teuchos::arcpFromArray(n_entries_per_row)); # endif // This functions assumes that the sparsity pattern sits on all @@ -150,7 +262,7 @@ namespace LinearAlgebra // that is already distributed. // now insert the indices - std::vector row_indices; + Teuchos::Array row_indices; for (const auto global_row : relevant_rows) { @@ -164,9 +276,7 @@ namespace LinearAlgebra sparsity_pattern.column_number(global_row, col); AssertIndexRange(global_row, row_space_map->getGlobalNumElements()); - graph->insertGlobalIndices(global_row, - row_length, - row_indices.data()); + graph->insertGlobalIndices(global_row, row_indices); } // Eventually, optimize the graph structure (sort indices, make memory -- 2.39.5