From 936ba3e20b2da6ddc96871c7bebc5c1c4f552c9c Mon Sep 17 00:00:00 2001 From: Sebastian Kinnewig Date: Sun, 18 Feb 2024 16:14:31 +0100 Subject: [PATCH] Add reinit() method to create a LA::TpetraWrappers::SparseMatrix from a dealii::SparseMatrix. --- .../lac/trilinos_tpetra_sparse_matrix.h | 28 +++++ .../trilinos_tpetra_sparse_matrix.templates.h | 102 ++++++++++++++++++ 2 files changed, 130 insertions(+) diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h index 543ceb8ec6..ad91ed457a 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h @@ -24,6 +24,8 @@ # include # include +# include +# include # include # include @@ -369,6 +371,32 @@ namespace LinearAlgebra const SparsityPatternType &sparsity_pattern, const MPI_Comm communicator = MPI_COMM_WORLD, const bool exchange_data = false); + + /** + * This function initializes the Trilinos matrix using the deal.II sparse + * matrix and the entries stored therein. It uses a threshold to copy only + * elements with modulus larger than the threshold (so zeros in the + * deal.II matrix can be filtered away). In contrast to the other reinit + * function with deal.II sparse matrix argument, this function takes a + * %parallel partitioning specified by the user instead of internally + * generating it. + * + * The optional parameter copy_values decides whether only the + * sparsity structure of the input matrix should be used or the matrix + * entries should be copied, too. + * + * This is a @ref GlossCollectiveOperation "collective operation" that needs to be called on all + * processors in order to avoid a dead lock. + */ + void + reinit(const IndexSet &row_parallel_partitioning, + const IndexSet &col_parallel_partitioning, + const dealii::SparseMatrix &dealii_sparse_matrix, + const MPI_Comm communicator = MPI_COMM_WORLD, + const double drop_tolerance = 1e-13, + const bool copy_values = true, + const dealii::SparsityPattern *use_this_sparsity = nullptr); + /** @} */ /** 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 650c1feb20..48b836e320 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h @@ -583,6 +583,108 @@ namespace LinearAlgebra + template + void + SparseMatrix::reinit( + const IndexSet &row_parallel_partitioning, + const IndexSet &col_parallel_partitioning, + const dealii::SparseMatrix &dealii_sparse_matrix, + const MPI_Comm communicator, + const double drop_tolerance, + const bool copy_values, + const dealii::SparsityPattern *use_this_sparsity) + { + dealii::types::signed_global_dof_index n_rows = dealii_sparse_matrix.m(); + AssertDimension(row_parallel_partitioning.size(), n_rows); + AssertDimension(col_parallel_partitioning.size(), + dealii_sparse_matrix.n()); + + const dealii::SparsityPattern &sparsity_pattern = + (use_this_sparsity != nullptr) ? + *use_this_sparsity : + dealii_sparse_matrix.get_sparsity_pattern(); + + if (matrix.is_null() || m() != n_rows || + n_nonzero_elements() != sparsity_pattern.n_nonzero_elements() || + copy_values) + if (use_this_sparsity == nullptr) + reinit(row_parallel_partitioning, + col_parallel_partitioning, + sparsity_pattern, + communicator, + false); + + // in case we do not copy values, we are done + if (copy_values == false) + return; + + // fill the values: go through all rows of the + // matrix, and then all columns. since the sparsity patterns of the + // input matrix and the specified sparsity pattern might be different, + // need to go through the row for both these sparsity structures + // simultaneously in order to really set the correct values. +# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0) + size_type maximum_row_length = matrix->getLocalMaxNumRowEntries(); +# else + size_type maximum_row_length = matrix->getNodeMaxNumRowEntries(); +# endif + std::vector row_indices(maximum_row_length); + std::vector values(maximum_row_length); + + for (dealii::types::signed_global_dof_index row = 0; row < n_rows; ++row) + // see if the row is locally stored on this processor + if (row_parallel_partitioning.is_element(row) == true) + { + dealii::SparsityPattern::iterator select_index = + sparsity_pattern.begin(row); + typename dealii::SparseMatrix::const_iterator it = + dealii_sparse_matrix.begin(row); + size_type col = 0; + if (sparsity_pattern.n_rows() == sparsity_pattern.n_cols()) + { + // optimized diagonal + AssertDimension(it->column(), row); + if (std::fabs(it->value()) > drop_tolerance) + { + values[col] = it->value(); + row_indices[col++] = it->column(); + } + ++select_index; + ++it; + } + + while (it != dealii_sparse_matrix.end(row) && + select_index != sparsity_pattern.end(row)) + { + while (select_index->column() < it->column() && + select_index != sparsity_pattern.end(row)) + ++select_index; + while (it->column() < select_index->column() && + it != dealii_sparse_matrix.end(row)) + ++it; + + if (it == dealii_sparse_matrix.end(row)) + break; + if (std::fabs(it->value()) > drop_tolerance) + { + values[col] = it->value(); + row_indices[col++] = it->column(); + } + ++select_index; + ++it; + } + set(row, + col, + reinterpret_cast(row_indices.data()), + values.data(), + false); + } + + compress(VectorOperation::insert); + } + + + // Information on the matrix template -- 2.39.5