Implement TrilinosWrappers::VectorBase::locally_owned_elements for vectors with non...

git-svn-id: https://svn.dealii.org/trunk@29638 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/lac/constraint_matrix.templates.h b/deal.II/include/deal.II/lac/constraint_matrix.templates.h
index 5cbed5d29f1e356005c216055dfcb363bf89ebc7..6e1fa5dba4eff544e0862e720a64bd4fa827cfcc 100644 (file)
--- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h
+++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h
@@ -745,6 +745,16 @@ namespace internal
             vec(i) = 0;
       }
 
+      // TODO: in general we should iterate over the constraints and not over all DoFs
+      // for performance reasons
+      template<typename Number>
+      void set_zero_parallel(const dealii::ConstraintMatrix &cm, parallel::distributed::Vector<Number> &vec, unsigned int shift = 0)
+      {
+        for (unsigned int i=0; i<vec.local_size(); ++i)
+          if (cm.is_constrained (shift + vec.local_range().first+i))
+            vec.local_element(i) = 0;
+      }
+
       template<class VEC>
       void set_zero_in_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, internal::bool2type<false>)
       {
@@ -794,21 +804,6 @@ namespace internal
       {
         set_zero_serial(cm, vec);
       }
-
-
-      template<class T>
-      void set_zero_all(const dealii::ConstraintMatrix &, dealii::parallel::distributed::Vector<T> &)
-      {
-        // should use the general template above, but requires that
-        // some member functions are added to parallel::distributed::Vector
-        Assert (false, ExcNotImplemented());
-      }
-
-      template<class T>
-      void set_zero_all(const dealii::ConstraintMatrix &, dealii::parallel::distributed::BlockVector<T> &)
-      {
-        Assert (false, ExcNotImplemented());
-      }
     }
   }
 }

diff --git a/deal.II/include/deal.II/lac/trilinos_vector_base.h b/deal.II/include/deal.II/lac/trilinos_vector_base.h
index 05d5f2411716271c932b32204cc68206ba133916..f8d2b168b387e682011263d9e2f530c26bedaf0a 100644 (file)
--- a/deal.II/include/deal.II/lac/trilinos_vector_base.h
+++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h
@@ -1248,11 +1248,27 @@ namespace TrilinosWrappers
   {
     IndexSet is (size());
 
-    // TODO: Trilinos does allow non-contiguous ranges for locally
-    // owned elements. if that is the case for the current
-    // vector, local_range() will throw an exception. Fix this.
-    const std::pair<unsigned int, unsigned int> x = local_range();
-    is.add_range (x.first, x.second);
+    // easy case: local range is contiguous
+    if (vector->Map().LinearMap())
+      {
+        const std::pair<unsigned int, unsigned int> x = local_range();
+        is.add_range (x.first, x.second);
+      }
+    else if (vector->Map().NumMyElements() > 0)
+      {
+        const unsigned int n_indices = vector->Map().NumMyElements();
+        int * vector_indices = vector->Map().MyGlobalElements();
+        unsigned int range_start = vector_indices[0];
+        unsigned int range_end = range_start+1;
+        for (unsigned int i=1; i<n_indices; ++i, ++range_end)
+          if (static_cast<unsigned int>(vector_indices[i]) != range_end)
+            {
+              is.add_range(range_start, range_end);
+              range_end = range_start = vector_indices[i];
+            }
+        is.add_range(range_start, range_end);
+        is.compress();
+      }
 
     return is;
   }

diff --git a/tests/mpi/trilinos_locally_owned_elements_01.cc b/tests/mpi/trilinos_locally_owned_elements_01.cc
new file mode 100644 (file)
index 0000000..7ed3dcf
--- /dev/null
+++ b/tests/mpi/trilinos_locally_owned_elements_01.cc
@@ -0,0 +1,88 @@
+//---------------------------------------------------------------------------
+//    $Id: trilinos_distribute_02.cc 29593 2013-05-25 12:29:14Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2009, 2010, 2012, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+// test TrilinosVector::locally_owned_elements
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <fstream>
+#include <sstream>
+
+
+
+template <int dim>
+void test()
+{
+  const unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  const unsigned int n_processes = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD);
+
+  // create non-contiguous index set
+  {
+    Assert(n_processes == 2, ExcNotImplemented());
+    IndexSet index (10);
+    for (unsigned int i=0; i<10; i+=2)
+      index.add_range(i+myid, i+myid+1);
+    index.compress();
+
+    TrilinosWrappers::MPI::Vector vec(index, MPI_COMM_WORLD);
+
+    IndexSet index2 = vec.locally_owned_elements();
+    Assert (index == index2, ExcInternalError());
+  }
+
+  // create contiguous index set
+  {
+    IndexSet index (10);
+    index.add_range(5*myid, 5*(myid+1));
+    index.compress();
+
+    TrilinosWrappers::MPI::Vector vec(index, MPI_COMM_WORLD);
+
+    IndexSet index2 = vec.locally_owned_elements();
+    Assert (index == index2, ExcInternalError());
+  }
+
+  if (myid == 0)
+    deallog << "OK" << std::endl;
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+
+
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      std::ofstream logfile(output_file_for_mpi("trilinos_locally_owned_elements_01").c_str());
+      deallog.attach(logfile);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+      test<2>();
+    }
+  else
+    {
+      deallog.depth_console(0);
+      test<2>();
+    }
+
+}

diff --git a/tests/mpi/trilinos_locally_owned_elements_01/ncpu_2/cmp/generic b/tests/mpi/trilinos_locally_owned_elements_01/ncpu_2/cmp/generic
new file mode 100644 (file)
index 0000000..be8d055
--- /dev/null
+++ b/tests/mpi/trilinos_locally_owned_elements_01/ncpu_2/cmp/generic
@@ -0,0 +1,2 @@
+
+DEAL:0::OK