Merge pull request #16672 from masterleinad/tpetra_is_non_negative

Tpetra: Implement Vector::is_non_negative

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;
+
       /** @} */
 
 

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<const Number> 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 <typename Number, typename MemorySpace>
+    bool
+    Vector<Number, MemorySpace>::is_non_negative() const
+    {
+      // get a representation of the vector and
+      // loop over all the elements
+      Teuchos::ArrayRCP<const Number> 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 <typename Number, typename MemorySpace>
     Number
     Vector<Number, MemorySpace>::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);

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<double>::is_non_zero
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::Vector<double> &v)
+{
+  // set only certain elements of the
+  // vector. they are all positive
+  std::vector<bool> 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<double> 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;
+    };
+}

tests/trilinos_tpetra/57.output

@@ -0,0 +1,2 @@
+
+DEAL::OK