From 4bd25fc569e079b48ee06b452878bae17de0735f Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 8 Feb 2024 23:39:34 -0500 Subject: [PATCH] Tpetra: Implement serial constructors, operator(), el(), and diag_element() for SparseMatrix --- .../lac/trilinos_tpetra_sparse_matrix.h | 95 ++++++++++++ .../trilinos_tpetra_sparse_matrix.templates.h | 140 ++++++++++++++++++ tests/trilinos_tpetra/01.cc | 100 +++++++++++++ tests/trilinos_tpetra/01.output | 2 + tests/trilinos_tpetra/02.cc | 106 +++++++++++++ tests/trilinos_tpetra/02.output | 2 + tests/trilinos_tpetra/03a.cc | 110 ++++++++++++++ tests/trilinos_tpetra/03a.output | 2 + tests/trilinos_tpetra/04.cc | 82 ++++++++++ tests/trilinos_tpetra/04.output | 2 + tests/trilinos_tpetra/05.cc | 96 ++++++++++++ tests/trilinos_tpetra/05.output | 3 + tests/trilinos_tpetra/09.cc | 105 +++++++++++++ tests/trilinos_tpetra/09.output | 2 + tests/trilinos_tpetra/CMakeLists.txt | 6 + 15 files changed, 853 insertions(+) create mode 100644 tests/trilinos_tpetra/01.cc create mode 100644 tests/trilinos_tpetra/01.output create mode 100644 tests/trilinos_tpetra/02.cc create mode 100644 tests/trilinos_tpetra/02.output create mode 100644 tests/trilinos_tpetra/03a.cc create mode 100644 tests/trilinos_tpetra/03a.output create mode 100644 tests/trilinos_tpetra/04.cc create mode 100644 tests/trilinos_tpetra/04.output create mode 100644 tests/trilinos_tpetra/05.cc create mode 100644 tests/trilinos_tpetra/05.output create mode 100644 tests/trilinos_tpetra/09.cc create mode 100644 tests/trilinos_tpetra/09.output create mode 100644 tests/trilinos_tpetra/CMakeLists.txt diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h index 772202b2e2..b145e7a220 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h @@ -152,6 +152,28 @@ namespace LinearAlgebra */ SparseMatrix(const SparsityPattern &sparsity_pattern); + /** + * Generate a matrix that is completely stored locally, having #m rows and + * #n columns. + * + * The number of columns entries per row is specified as the maximum + * number of entries argument. + */ + SparseMatrix(const size_type m, + const size_type n, + const unsigned int n_max_entries_per_row); + + /** + * Generate a matrix that is completely stored locally, having #m rows and + * #n columns. + * + * The vector n_entries_per_row specifies the number of entries + * in each row. + */ + SparseMatrix(const size_type m, + const size_type n, + const std::vector &n_entries_per_row); + /** * Move constructor. Create a new sparse matrix by stealing the internal * data of the `other` object. @@ -587,7 +609,50 @@ namespace LinearAlgebra const bool elide_zero_values = false); /** @} */ + /** + * @name Entry Access + */ + /** @{ */ + /** + * Return the value of the entry (i,j). This may be an expensive + * operation and you should always take care where to call this function. + * As in the deal.II sparse matrix class, we throw an exception if the + * respective entry doesn't exist in the sparsity pattern of this class, + * which is requested from Trilinos. Moreover, an exception will be thrown + * when the requested element is not saved on the calling process. + */ + Number + operator()(const size_type i, const size_type j) const; + + /** + * Return the value of the matrix entry (i,j). If this entry does + * not exist in the sparsity pattern, then zero is returned. While this + * may be convenient in some cases, note that it is simple to write + * algorithms that are slow compared to an optimal solution, since the + * sparsity of the matrix is not used. On the other hand, if you want to + * be sure the entry exists, you should use operator() instead. + * + * The lack of error checking in this function can also yield surprising + * results if you have a parallel matrix. In that case, just because you + * get a zero result from this function does not mean that either the + * entry does not exist in the sparsity pattern or that it does but has a + * value of zero. Rather, it could also be that it simply isn't stored on + * the current processor; in that case, it may be stored on a different + * processor, and possibly so with a nonzero value. + */ + Number + el(const size_type i, const size_type j) const; + + /** + * Return the main diagonal element in the ith row. This function + * throws an error if the matrix is not quadratic and it also throws an + * error if (i,i) is not element of the local matrix. + */ + Number + diag_element(const size_type i) const; + + /** @} */ /** * @name Multiplications */ @@ -802,9 +867,39 @@ namespace LinearAlgebra "Are you trying to put the product of the " "matrix with a vector into a vector that has " "ghost elements?"); + + /** + * Exception + */ + DeclException2(ExcInvalidIndex, + size_type, + size_type, + << "The entry with index <" << arg1 << ',' << arg2 + << "> does not exist."); + + /** + * Exception + */ + DeclException4(ExcAccessToNonLocalElement, + size_type, + size_type, + size_type, + size_type, + << "You tried to access element (" << arg1 << '/' << arg2 + << ')' + << " of a distributed matrix, but only rows in range [" + << arg3 << ',' << arg4 + << "] are stored locally and can be accessed."); + /** @} */ private: + /** + * Helper function for operator() and el(). + */ + Number + element(const size_type i, const size_type j, const bool no_error) const; + /** * Pointer to the user-supplied Tpetra Trilinos mapping of the matrix * columns that assigns parts of the matrix to the individual processes. 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 6e9ae879b5..8273703b26 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h @@ -230,6 +230,49 @@ namespace LinearAlgebra + template + SparseMatrix::SparseMatrix( + const size_type m, + const size_type n, + const unsigned int n_max_entries_per_row) + : column_space_map(Utilities::Trilinos::internal::make_rcp( + n, + 0, + Utilities::Trilinos::tpetra_comm_self())) + , matrix(Utilities::Trilinos::internal::make_rcp( + Utilities::Trilinos::internal::make_rcp( + m, + 0, + Utilities::Trilinos::tpetra_comm_self()), + column_space_map, + n_max_entries_per_row)) + , compressed(false) + {} + + + + template + SparseMatrix::SparseMatrix( + const size_type m, + const size_type n, + const std::vector &n_entries_per_row) + : column_space_map(Utilities::Trilinos::internal::make_rcp( + n, + 0, + Utilities::Trilinos::tpetra_comm_self())) + , compressed(false) + { + std::vector entries_per_row_size_type(n_entries_per_row.begin(), + n_entries_per_row.end()); + matrix = Utilities::Trilinos::internal::make_rcp( + Utilities::Trilinos::internal::make_rcp( + m, 0, Utilities::Trilinos::tpetra_comm_self()), + column_space_map, + Teuchos::ArrayView{entries_per_row_size_type}); + } + + + template SparseMatrix::SparseMatrix( SparseMatrix &&other) noexcept @@ -826,6 +869,8 @@ namespace LinearAlgebra } } + + template void SparseMatrix::resume_fill() @@ -837,6 +882,101 @@ namespace LinearAlgebra } } + + template + Number + SparseMatrix::element(const size_type i, + const size_type j, + const bool no_error) const + { + // Extract local indices in the matrix. + const int trilinos_i = matrix->getRowMap()->getLocalElement(i); + const int trilinos_j = matrix->getColMap()->getLocalElement(j); + TrilinosScalar value = 0.; + + if (trilinos_i == Teuchos::OrdinalTraits::invalid() || + trilinos_j == Teuchos::OrdinalTraits::invalid()) + { + if (no_error) + return {}; + Assert(false, + ExcAccessToNonLocalElement( + i, j, local_range().first, local_range().second - 1)); + } + else + { + Assert(matrix->isFillComplete(), ExcMatrixNotCompressed()); + + // Prepare pointers for extraction of a view of the row. + size_t nnz_present = matrix->getNumEntriesInLocalRow(trilinos_i); + typename MatrixType::nonconst_local_inds_host_view_type col_indices( + "indices", nnz_present); + typename MatrixType::nonconst_values_host_view_type values( + "values", nnz_present); + + matrix->getLocalRowCopy(trilinos_i, col_indices, values, nnz_present); + + // Search the index where we look for the value, and then finally get + // it. + int local_col_index = 0; + for (; local_col_index < static_cast(nnz_present); + ++local_col_index) + { + if (col_indices(local_col_index) == trilinos_j) + break; + } + + if (local_col_index == static_cast(nnz_present)) + { + if (no_error) + return {}; + Assert(false, ExcInvalidIndex(i, j)); + } + else + value = values(local_col_index); + } + + return value; + } + + + + template + Number + SparseMatrix::operator()(const size_type i, + const size_type j) const + { + return element(i, j, /* no_error */ false); + } + + + + template + Number + SparseMatrix::el(const size_type i, + const size_type j) const + { + return element(i, j, /* no_error */ true); + } + + + + template + Number + SparseMatrix::diag_element(const size_type i) const + { + Assert(m() == n(), ExcNotQuadratic()); + +# ifdef DEBUG + // use operator() in debug mode because it checks if this is a valid + // element (in parallel) + return operator()(i, i); +# else + // Trilinos doesn't seem to have a more efficient way to access the + // diagonal than by just using the standard el(i,j) function. + return el(i, i); +# endif + } } // namespace TpetraWrappers } // namespace LinearAlgebra diff --git a/tests/trilinos_tpetra/01.cc b/tests/trilinos_tpetra/01.cc new file mode 100644 index 0000000000..e39faabba7 --- /dev/null +++ b/tests/trilinos_tpetra/01.cc @@ -0,0 +1,100 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check setting elements in a trilinos matrix using +// LinearAlgebra::TpetraWrappers::SparseMatrix::set() + +#include + +#include + +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + // first set a few entries + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + m.set(i, j, i * j * .5 + .5); + + m.compress(VectorOperation::insert); + + // then make sure we retrieve the same ones + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + { + AssertThrow(m(i, j) == i * j * .5 + .5, ExcInternalError()); + AssertThrow(m.el(i, j) == i * j * .5 + .5, ExcInternalError()); + } + else + { + AssertThrow(m.el(i, j) == 0, ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix m(5U, 5U, 3U); + test(m); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/01.output b/tests/trilinos_tpetra/01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos_tpetra/01.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/trilinos_tpetra/02.cc b/tests/trilinos_tpetra/02.cc new file mode 100644 index 0000000000..e620df8f69 --- /dev/null +++ b/tests/trilinos_tpetra/02.cc @@ -0,0 +1,106 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check setting elements in a petsc matrix using +// LinearAlgebra::TpetraWrappers::SparseMatrix::add() + +#include + +#include + +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + // first set a few entries + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + m.set(i, j, 0.); + + // then add values to these entries + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + m.add(i, j, i * j * .5 + .5); + + m.compress(VectorOperation::add); + + // then make sure we retrieve the same ones + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + { + AssertThrow(m(i, j) == i * j * .5 + .5, ExcInternalError()); + AssertThrow(m.el(i, j) == i * j * .5 + .5, ExcInternalError()); + } + else + { + AssertThrow(m.el(i, j) == 0, ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix m(5U, 5U, 3U); + test(m); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/02.output b/tests/trilinos_tpetra/02.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos_tpetra/02.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/trilinos_tpetra/03a.cc b/tests/trilinos_tpetra/03a.cc new file mode 100644 index 0000000000..3c63425542 --- /dev/null +++ b/tests/trilinos_tpetra/03a.cc @@ -0,0 +1,110 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check setting elements in a petsc matrix using set() and add() +// intermixed. this poses PETSc some problems, since one has to flush some +// buffer in between these two types of operations +// +// in contrast to trilinos_03, we set and add the same elements here twice, to +// get double the original value + +#include + +#include + +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + // first set a few entries + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + m.set(i, j, i * j * .5 + .5); + + // then add the same elements again + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + m.add(i, j, i * j * .5 + .5); + + m.compress(VectorOperation::add); + + // then make sure we retrieve the same ones + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + { + AssertThrow(m(i, j) == 2 * (i * j * .5 + .5), ExcInternalError()); + AssertThrow(m.el(i, j) == 2 * (i * j * .5 + .5), ExcInternalError()); + } + else + { + AssertThrow(m.el(i, j) == 0, ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix m(5U, 5U, 3U); + test(m); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/03a.output b/tests/trilinos_tpetra/03a.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos_tpetra/03a.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/trilinos_tpetra/04.cc b/tests/trilinos_tpetra/04.cc new file mode 100644 index 0000000000..86a1942383 --- /dev/null +++ b/tests/trilinos_tpetra/04.cc @@ -0,0 +1,82 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check querying matrix sizes + +#include + +#include + +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + AssertThrow(m.m() == 5, ExcInternalError()); + AssertThrow(m.n() == 5, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix m(5U, 5U, 3U); + test(m); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/04.output b/tests/trilinos_tpetra/04.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos_tpetra/04.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/trilinos_tpetra/05.cc b/tests/trilinos_tpetra/05.cc new file mode 100644 index 0000000000..7561408b65 --- /dev/null +++ b/tests/trilinos_tpetra/05.cc @@ -0,0 +1,96 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check querying the number of nonzero elements in +// LinearAlgebra::TpetraWrappers::SparseMatrix + +#include + +#include + +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + // first set a few entries. count how many + // entries we have + unsigned int counter = 0; + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + { + m.set(i, j, i * j * .5 + .5); + ++counter; + } + + m.compress(VectorOperation::insert); + + deallog << m.n_nonzero_elements() << std::endl; + Assert(m.n_nonzero_elements() == counter, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix m(5U, 5U, 3U); + test(m); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/05.output b/tests/trilinos_tpetra/05.output new file mode 100644 index 0000000000..169db65e4f --- /dev/null +++ b/tests/trilinos_tpetra/05.output @@ -0,0 +1,3 @@ + +DEAL::8 +DEAL::OK diff --git a/tests/trilinos_tpetra/09.cc b/tests/trilinos_tpetra/09.cc new file mode 100644 index 0000000000..4016f89e05 --- /dev/null +++ b/tests/trilinos_tpetra/09.cc @@ -0,0 +1,105 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check LinearAlgebra::TpetraWrappers::SparseMatrix::operator *= + +#include + +#include + +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + // first set a few entries + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + m.set(i, j, i * j * .5 + .5); + + m.compress(VectorOperation::insert); + + // then multiply everything by 1.25 and + // make sure we retrieve the values we + // expect + m *= 1.25; + + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.m(); ++j) + if ((i + 2 * j + 1) % 3 == 0) + { + AssertThrow(m(i, j) == (i * j * .5 + .5) * 1.25, ExcInternalError()); + AssertThrow(m.el(i, j) == (i * j * .5 + .5) * 1.25, + ExcInternalError()); + } + else + { + AssertThrow(m.el(i, j) == 0, ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix m(5U, 5U, 3U); + test(m); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/09.output b/tests/trilinos_tpetra/09.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos_tpetra/09.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/trilinos_tpetra/CMakeLists.txt b/tests/trilinos_tpetra/CMakeLists.txt new file mode 100644 index 0000000000..1a131cbb32 --- /dev/null +++ b/tests/trilinos_tpetra/CMakeLists.txt @@ -0,0 +1,6 @@ +cmake_minimum_required(VERSION 3.13.4) +include(../scripts/setup_testsubproject.cmake) +project(testsuite CXX) +if(DEAL_II_WITH_TRILINOS) + deal_ii_pickup_tests() +endif() -- 2.39.5