From: Daniel Arndt Date: Tue, 6 Feb 2024 17:58:48 +0000 (-0500) Subject: Tpetra: Enable lac/linear_operator_15 X-Git-Tag: relicensing~67^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6d0306c355c2ea4a63b601987c59425b77acac88;p=dealii.git Tpetra: Enable lac/linear_operator_15 --- diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h index bbaa1e3ba8..772202b2e2 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h @@ -446,6 +446,146 @@ namespace LinearAlgebra const TrilinosScalar *values, const bool elide_zero_values = true, const bool col_indices_are_sorted = false); + + /** + * Set the element (i,j) to @p value. + * + * This function is able to insert new elements into the matrix as long as + * compress() has not been called, so the sparsity pattern will be + * extended. When compress() is called for the first time (or in case the + * matrix is initialized from a sparsity pattern), no new elements can be + * added and an insertion of elements at positions which have not been + * initialized will throw an exception. + * + * For the case that the matrix is constructed without a sparsity pattern + * and new matrix entries are added on demand, please note the following + * behavior imposed by the underlying Tpetra::CrsMatrix data structure: + * If the same matrix entry is inserted more than once, the matrix entries + * will be added upon calling compress() (since Tpetra does not track + * values to the same entry before the final compress() is called), even + * if VectorOperation::insert is specified as argument to compress(). In + * the case you cannot make sure that matrix entries are only set once, + * initialize the matrix with a sparsity pattern to fix the matrix + * structure before inserting elements. + */ + void + set(const size_type i, const size_type j, const Number value); + + /** + * Set all elements given in a FullMatrix into the sparse matrix + * locations given by indices. In other words, this function + * writes the elements in full_matrix into the calling matrix, + * using the local-to-global indexing specified by indices for + * both the rows and the columns of the matrix. This function assumes a + * quadratic sparse matrix and a quadratic full_matrix, the usual + * situation in FE calculations. + * + * This function is able to insert new elements into the matrix as long as + * compress() has not been called, so the sparsity pattern will be + * extended. After compress() has been called for the first time or the + * matrix has been initialized from a sparsity pattern, extending the + * sparsity pattern is no longer possible and an insertion of elements at + * positions which have not been initialized will throw an exception. + * + * The optional parameter elide_zero_values can be used to + * specify whether zero values should be inserted anyway or they should be + * filtered away. The default value is false, i.e., even zero + * values are inserted/replaced. + * + * For the case that the matrix is constructed without a sparsity pattern + * and new matrix entries are added on demand, please note the following + * behavior imposed by the underlying Tpetra::CrsMatrix data structure: + * If the same matrix entry is inserted more than once, the matrix entries + * will be added upon calling compress() (since Epetra does not track + * values to the same entry before the final compress() is called), even + * if VectorOperation::insert is specified as argument to compress(). In + * the case you cannot make sure that matrix entries are only set once, + * initialize the matrix with a sparsity pattern to fix the matrix + * structure before inserting elements. + */ + void + set(const std::vector &indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = false); + + /** + * Same function as before, but now including the possibility to use + * rectangular full_matrices and different local-to-global indexing on + * rows and columns, respectively. + */ + void + set(const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = false); + + /** + * Set several elements in the specified row of the matrix with column + * indices as given by col_indices to the respective value. + * + * This function is able to insert new elements into the matrix as long as + * compress() has not been called, so the sparsity pattern will be + * extended. After compress() has been called for the first time or the + * matrix has been initialized from a sparsity pattern, extending the + * sparsity pattern is no longer possible and an insertion of elements at + * positions which have not been initialized will throw an exception. + * + * The optional parameter elide_zero_values can be used to + * specify whether zero values should be inserted anyway or they should be + * filtered away. The default value is false, i.e., even zero + * values are inserted/replaced. + * + * For the case that the matrix is constructed without a sparsity pattern + * and new matrix entries are added on demand, please note the following + * behavior imposed by the underlying Tpetra::CrsMatrix data structure: + * If the same matrix entry is inserted more than once, the matrix entries + * will be added upon calling compress() (since Epetra does not track + * values to the same entry before the final compress() is called), even + * if VectorOperation::insert is specified as argument to compress(). In + * the case you cannot make sure that matrix entries are only set once, + * initialize the matrix with a sparsity pattern to fix the matrix + * structure before inserting elements. + */ + void + set(const size_type row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values = false); + + /** + * Set several elements to values given by values in a given row + * in columns given by col_indices into the sparse matrix. + * + * This function is able to insert new elements into the matrix as long as + * compress() has not been called, so the sparsity pattern will be + * extended. After compress() has been called for the first time or the + * matrix has been initialized from a sparsity pattern, extending the + * sparsity pattern is no longer possible and an insertion of elements at + * positions which have not been initialized will throw an exception. + * + * The optional parameter elide_zero_values can be used to + * specify whether zero values should be inserted anyway or they should be + * filtered away. The default value is false, i.e., even zero + * values are inserted/replaced. + * + * For the case that the matrix is constructed without a sparsity pattern + * and new matrix entries are added on demand, please note the following + * behavior imposed by the underlying Tpetra::CrsMatrix data structure: + * If the same matrix entry is inserted more than once, the matrix entries + * will be added upon calling compress() (since Epetra does not track + * values to the same entry before the final compress() is called), even + * if VectorOperation::insert is specified as argument to compress(). In + * the case you cannot make sure that matrix entries are only set once, + * initialize the matrix with a sparsity pattern to fix the matrix + * structure before inserting elements. + */ + void + set(const size_type row, + const size_type n_cols, + const size_type *col_indices, + const Number *values, + const bool elide_zero_values = false); + /** @} */ /** 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 9cb4244606..6e9ae879b5 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h @@ -21,6 +21,7 @@ #ifdef DEAL_II_TRILINOS_WITH_TPETRA # include +# include # include # include @@ -516,6 +517,118 @@ namespace LinearAlgebra + template + void + SparseMatrix::set(const size_type row, + const size_type n_cols, + const size_type *col_indices, + const Number *values, + const bool elide_zero_values) + { + AssertIndexRange(row, this->m()); + const types::signed_global_dof_index *col_index_ptr; + const Number *col_value_ptr; + const types::signed_global_dof_index trilinos_row = row; + types::signed_global_dof_index n_columns; + + boost::container::small_vector local_value_array( + elide_zero_values ? n_cols : 0); + boost::container::small_vector + local_index_array(elide_zero_values ? n_cols : 0); + + // If we don't elide zeros, the pointers are already available... need to + // cast to non-const pointers as that is the format taken by Trilinos (but + // we will not modify const data) + if (elide_zero_values == false) + { + col_index_ptr = + reinterpret_cast( + col_indices); + col_value_ptr = values; + n_columns = n_cols; + } + else + { + // Otherwise, extract nonzero values in each row and get the + // respective indices. + col_index_ptr = local_index_array.data(); + col_value_ptr = local_value_array.data(); + + n_columns = 0; + for (size_type j = 0; j < n_cols; ++j) + { + const double value = values[j]; + AssertIsFinite(value); + if (value != 0) + { + local_index_array[n_columns] = col_indices[j]; + local_value_array[n_columns] = value; + ++n_columns; + } + } + + AssertIndexRange(n_columns, n_cols + 1); + } + + // We distinguish between two cases: the first one is when the matrix is + // not filled (i.e., it is possible to add new elements to the sparsity + // pattern), and the second one is when the pattern is already fixed. In + // the former case, we add the possibility to insert new values, and in + // the second we just replace data. + if (compressed || matrix->isFillComplete()) + { + matrix->resumeFill(); + compressed = false; + } + + if (!matrix->isStaticGraph()) + matrix->insertGlobalValues(trilinos_row, + n_columns, + col_value_ptr, + col_index_ptr); + else + matrix->replaceGlobalValues(trilinos_row, + n_columns, + col_value_ptr, + col_index_ptr); + } + + + + template + inline void + SparseMatrix::set(const size_type i, + const size_type j, + const Number value) + { + AssertIsFinite(value); + + set(i, 1, &j, &value, false); + } + + + + template + inline void + SparseMatrix::set( + const std::vector &indices, + const FullMatrix &values, + const bool elide_zero_values) + { + Assert(indices.size() == values.m(), + ExcDimensionMismatch(indices.size(), values.m())); + Assert(values.m() == values.n(), ExcNotQuadratic()); + + for (size_type i = 0; i < indices.size(); ++i) + set(indices[i], + indices.size(), + indices.data(), + &values(i, 0), + elide_zero_values); + } + + + template void SparseMatrix::add( diff --git a/include/deal.II/lac/trilinos_tpetra_vector.h b/include/deal.II/lac/trilinos_tpetra_vector.h index 80a140b1c7..9ac1bd8cc2 100644 --- a/include/deal.II/lac/trilinos_tpetra_vector.h +++ b/include/deal.II/lac/trilinos_tpetra_vector.h @@ -335,7 +335,7 @@ namespace LinearAlgebra */ void reinit(const IndexSet ¶llel_partitioner, - const MPI_Comm communicator, + const MPI_Comm communicator = MPI_COMM_WORLD, const bool omit_zeroing_entries = false); /** diff --git a/source/lac/trilinos_tpetra_sparse_matrix.cc b/source/lac/trilinos_tpetra_sparse_matrix.cc index 6e5ce3fcc3..63f970f3ed 100644 --- a/source/lac/trilinos_tpetra_sparse_matrix.cc +++ b/source/lac/trilinos_tpetra_sparse_matrix.cc @@ -37,6 +37,9 @@ namespace LinearAlgebra const MPI_Comm communicator, const bool exchange_data); + template void + SparseMatrix::reinit(const dealii::DynamicSparsityPattern &); + } // namespace TpetraWrappers } // namespace LinearAlgebra # endif // DOXYGEN diff --git a/tests/lac/linear_operator_15.cc b/tests/lac/linear_operator_15.cc index f51d4553af..4a1efe3f55 100644 --- a/tests/lac/linear_operator_15.cc +++ b/tests/lac/linear_operator_15.cc @@ -22,6 +22,10 @@ #include #include #include +#ifdef DEAL_II_TRILINOS_WITH_TPETRA +# include +# include +#endif #include #include "../tests.h" @@ -91,6 +95,39 @@ main(int argc, char *argv[]) u = lo_A_plus_id_2 * f; } +#ifdef DEAL_II_TRILINOS_WITH_TPETRA + { + DynamicSparsityPattern csp(dim, dim); + testproblem.five_point_structure(csp); + + LinearAlgebra::TpetraWrappers::SparseMatrix A; + A.reinit(csp); + testproblem.five_point(A); + + LinearAlgebra::TpetraWrappers::Vector f; + f.reinit(complete_index_set(dim)); + LinearAlgebra::TpetraWrappers::Vector u; + u.reinit(complete_index_set(dim)); + + A.compress(VectorOperation::insert); + f.compress(VectorOperation::insert); + u.compress(VectorOperation::insert); + + const auto lo_A = + linear_operator>(A); + const auto lo_id_1 = + identity_operator>( + lo_A.reinit_range_vector); + const auto lo_id_2 = + identity_operator>(lo_A); + const auto lo_A_plus_id_1 = lo_A + lo_id_1; // Not a good idea. See below. + const auto lo_A_plus_id_2 = lo_A + lo_id_2; + + // u = lo_A_plus_id_1 * f; // This is not allowed -- different payloads + u = lo_A_plus_id_2 * f; + } +#endif + deallog << "OK" << std::endl; return 0;