From: Daniel Arndt Date: Mon, 19 Feb 2024 18:25:09 +0000 (-0500) Subject: Tpetra: Implement Vector::is_non_negative X-Git-Tag: relicensing~19^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=64b48bdc5790b220eae41e9d58d3c90525fcc30c;p=dealii.git Tpetra: Implement Vector::is_non_negative --- diff --git a/include/deal.II/lac/trilinos_tpetra_vector.h b/include/deal.II/lac/trilinos_tpetra_vector.h index e8bf4ca7d5..c2a90e7210 100644 --- a/include/deal.II/lac/trilinos_tpetra_vector.h +++ b/include/deal.II/lac/trilinos_tpetra_vector.h @@ -653,6 +653,14 @@ namespace LinearAlgebra bool all_zero() const; + /** + * Return @p true if the vector has no negative entries, i.e. all entries + * are zero or positive. This function is used, for example, to check + * whether refinement indicators are really all positive (or zero). + */ + bool + is_non_negative() const; + /** @} */ diff --git a/include/deal.II/lac/trilinos_tpetra_vector.templates.h b/include/deal.II/lac/trilinos_tpetra_vector.templates.h index 1f2aca18b0..04fd620478 100644 --- a/include/deal.II/lac/trilinos_tpetra_vector.templates.h +++ b/include/deal.II/lac/trilinos_tpetra_vector.templates.h @@ -808,18 +808,16 @@ namespace LinearAlgebra { // get a representation of the vector and // loop over all the elements - Number *start_ptr = vector->getDataNonConst().get(); - const Number *ptr = start_ptr, - *eptr = start_ptr + vector->getLocalLength(); - unsigned int flag = 0; - while (ptr != eptr) + Teuchos::ArrayRCP data = vector->getData(); + const size_type n_elements = vector->getLocalLength(); + unsigned int flag = 0; + for (size_type i = 0; i < n_elements; ++i) { - if (*ptr != Number(0)) + if (data[i] != Number(0)) { flag = 1; break; } - ++ptr; } // Check that the vector is zero on _all_ processors. @@ -833,6 +831,34 @@ namespace LinearAlgebra + template + bool + Vector::is_non_negative() const + { + // get a representation of the vector and + // loop over all the elements + Teuchos::ArrayRCP data = vector->getData(); + const size_type n_elements = vector->getLocalLength(); + unsigned int flag = 0; + for (size_type i = 0; i < n_elements; ++i) + { + if (data[i] < Number(0)) + { + flag = 1; + break; + } + } + + // Check that the vector is non-negative on _all_ processors. + unsigned int num_negative = + Utilities::MPI::sum(flag, + Utilities::Trilinos::teuchos_comm_to_mpi_comm( + vector->getMap()->getComm())); + return num_negative == 0; + } + + + template Number Vector::mean_value() const @@ -1066,10 +1092,6 @@ namespace LinearAlgebra AssertThrow(out.fail() == false, ExcIO()); boost::io::ios_flags_saver restore_flags(out); - // Get a representation of the vector and loop over all - // the elements - const auto val = vector->get1dView(); - out.precision(precision); if (scientific) out.setf(std::ios::scientific, std::ios::floatfield); diff --git a/tests/trilinos_tpetra/57.cc b/tests/trilinos_tpetra/57.cc new file mode 100644 index 0000000000..5123bf5780 --- /dev/null +++ b/tests/trilinos_tpetra/57.cc @@ -0,0 +1,101 @@ +// --------------------------------------------------------------------- +// +// 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::Vector::is_non_zero + +#include + +#include + +#include +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::Vector &v) +{ + // set only certain elements of the + // vector. they are all positive + std::vector pattern(v.size(), false); + for (unsigned int i = 0; i < v.size(); i += 1 + i) + { + v(i) += i; + pattern[i] = true; + } + + v.compress(VectorOperation::add); + + // check that the vector is really + // non-negative + AssertThrow(v.is_non_negative() == true, ExcInternalError()); + + // then set a single element to a negative + // value and check again + v(v.size() / 2) = -1; + AssertThrow(v.is_non_negative() == false, 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::Vector v; + v.reinit(complete_index_set(100), MPI_COMM_WORLD); + test(v); + } + } + 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/57.output b/tests/trilinos_tpetra/57.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos_tpetra/57.output @@ -0,0 +1,2 @@ + +DEAL::OK