From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Mon, 15 Aug 2016 15:30:09 +0000 (+0200)
Subject: Fix TrilinosWrappers::locally_owned_elements
X-Git-Tag: v8.5.0-rc1~667^2~5
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f94596e7f62faf15afcfa10972a45de7badca4ce;p=dealii.git

Fix TrilinosWrappers::locally_owned_elements
---

diff --git a/include/deal.II/lac/trilinos_vector.h b/include/deal.II/lac/trilinos_vector.h
index ef4cbb8783..8e57c44931 100644
--- a/include/deal.II/lac/trilinos_vector.h
+++ b/include/deal.II/lac/trilinos_vector.h
@@ -661,6 +661,7 @@ namespace TrilinosWrappers
     {
       *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
                      v);
+      owned_elements = parallel_partitioner;
     }
 
 
@@ -671,9 +672,7 @@ namespace TrilinosWrappers
                          const dealii::Vector<number> &v)
     {
       if (vector.get() == 0 || vector->Map().SameAs(parallel_partitioner) == false)
-        vector.reset (new Epetra_FEVector(parallel_partitioner));
-
-      has_ghosts = vector->Map().UniqueGIDs()==false;
+        reinit(parallel_partitioner);
 
       const int size = parallel_partitioner.NumMyElements();
 
@@ -681,6 +680,7 @@ namespace TrilinosWrappers
       // that a direct access is not possible.
       for (int i=0; i<size; ++i)
         (*vector)[0][i] = v(gid(parallel_partitioner,i));
+
     }
 
 
@@ -966,6 +966,7 @@ namespace TrilinosWrappers
         Epetra_LocalMap map ((TrilinosWrappers::types::int_type)v.size(), 0,
                              Utilities::Trilinos::comm_self());
         vector.reset (new Epetra_FEVector(map));
+        owned_elements = v.locally_owned_elements();
       }
 
     const Epetra_Map &map = vector_partitioner();
diff --git a/include/deal.II/lac/trilinos_vector_base.h b/include/deal.II/lac/trilinos_vector_base.h
index 1e71b479c0..a37bbe9613 100644
--- a/include/deal.II/lac/trilinos_vector_base.h
+++ b/include/deal.II/lac/trilinos_vector_base.h
@@ -937,6 +937,11 @@ namespace TrilinosWrappers
      */
     std_cxx11::shared_ptr<Epetra_MultiVector> nonlocal_vector;
 
+    /**
+     * An IndexSet storing the indices this vector owns exclusively.
+     */
+    IndexSet owned_elements;
+
     /**
      * Make the reference class a friend.
      */
@@ -1082,27 +1087,7 @@ namespace TrilinosWrappers
   IndexSet
   VectorBase::locally_owned_elements() const
   {
-    IndexSet is (size());
-
-    // easy case: local range is contiguous
-    if (vector->Map().LinearMap())
-      {
-        const std::pair<size_type, size_type> x = local_range();
-        is.add_range (x.first, x.second);
-      }
-    else if (vector->Map().NumMyElements() > 0)
-      {
-        const size_type n_indices = vector->Map().NumMyElements();
-#ifndef DEAL_II_WITH_64BIT_INDICES
-        unsigned int *vector_indices = (unsigned int *)vector->Map().MyGlobalElements();
-#else
-        size_type *vector_indices = (size_type *)vector->Map().MyGlobalElements64();
-#endif
-        is.add_indices(vector_indices, vector_indices+n_indices);
-        is.compress();
-      }
-
-    return is;
+    return owned_elements;
   }
 
 
@@ -1218,6 +1203,8 @@ namespace TrilinosWrappers
 
     if (v.nonlocal_vector.get() != 0)
       nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
+
+    owned_elements = v.locally_owned_elements();
   }
 
 
diff --git a/source/lac/trilinos_vector.cc b/source/lac/trilinos_vector.cc
index 9bdddcb538..7c3f6868e3 100644
--- a/source/lac/trilinos_vector.cc
+++ b/source/lac/trilinos_vector.cc
@@ -122,6 +122,7 @@ namespace TrilinosWrappers
       last_action = Zero;
       vector.reset (new Epetra_FEVector(*v.vector));
       has_ghosts = v.has_ghosts;
+      owned_elements = v.locally_owned_elements();
     }
 
 
@@ -132,6 +133,7 @@ namespace TrilinosWrappers
       // initialize a minimal, valid object and swap
       last_action = Zero;
       vector.reset(new Epetra_FEVector(Epetra_Map(0,0,0,Utilities::Trilinos::comm_self())));
+      owned_elements.clear();
 
       swap(v);
     }
@@ -186,9 +188,7 @@ namespace TrilinosWrappers
       :
       VectorBase()
     {
-      IndexSet parallel_partitioning = local;
-      parallel_partitioning.add_indices(ghost);
-      reinit(parallel_partitioning, communicator);
+      reinit(local, ghost, communicator, false);
     }
 
 
@@ -207,6 +207,30 @@ namespace TrilinosWrappers
       vector.reset (new Epetra_FEVector(input_map));
 
       has_ghosts = vector->Map().UniqueGIDs()==false;
+
+      {
+        owned_elements.clear();
+        owned_elements.set_size(size());
+
+        // easy case: local range is contiguous
+        if (vector->Map().LinearMap())
+          {
+            const std::pair<size_type, size_type> x = local_range();
+            owned_elements.add_range (x.first, x.second);
+          }
+        else if (vector->Map().NumMyElements() > 0)
+          {
+            const size_type n_indices = vector->Map().NumMyElements();
+#ifndef DEAL_II_WITH_64BIT_INDICES
+            unsigned int *vector_indices = (unsigned int *)vector->Map().MyGlobalElements();
+#else
+            size_type *vector_indices = (size_type *)vector->Map().MyGlobalElements64();
+#endif
+            owned_elements.add_indices(vector_indices, vector_indices+n_indices);
+            owned_elements.compress();
+          }
+      }
+
       last_action = Zero;
     }
 
@@ -215,13 +239,20 @@ namespace TrilinosWrappers
     void
     Vector::reinit (const IndexSet &parallel_partitioner,
                     const MPI_Comm &communicator,
-                    const bool      omit_zeroing_entries)
+                    const bool      /*omit_zeroing_entries*/)
     {
       nonlocal_vector.reset();
 
       Epetra_Map map = parallel_partitioner.make_trilinos_map (communicator,
                                                                true);
-      reinit (map, omit_zeroing_entries);
+
+      vector.reset (new Epetra_FEVector(map));
+
+      has_ghosts = vector->Map().UniqueGIDs()==false;
+
+      owned_elements = parallel_partitioner;
+
+      last_action = Zero;
     }
 
 
@@ -255,6 +286,7 @@ namespace TrilinosWrappers
               vector.reset (new Epetra_FEVector(v.vector->Map()));
               has_ghosts = v.has_ghosts;
               last_action = Zero;
+              owned_elements = v.locally_owned_elements();
             }
           else if (omit_zeroing_entries == false)
             {
@@ -303,6 +335,7 @@ namespace TrilinosWrappers
                     const bool         import_data)
     {
       nonlocal_vector.reset();
+      owned_elements.clear();
 
       // In case we do not allow to have different maps, this call means that
       // we have to reset the vector. So clear the vector, initialize our map
@@ -316,6 +349,7 @@ namespace TrilinosWrappers
       for (size_type block=0; block<v.n_blocks(); ++block)
         n_elements += v.block(block).local_size();
       std::vector<TrilinosWrappers::types::int_type> global_ids (n_elements, -1);
+      size_type max_size = 0;
       for (size_type block=0; block<v.n_blocks(); ++block)
         {
           TrilinosWrappers::types::int_type *glob_elements =
@@ -323,7 +357,9 @@ namespace TrilinosWrappers
           for (size_type i=0; i<v.block(block).local_size(); ++i)
             global_ids[added_elements++] = glob_elements[i] + block_offset;
           block_offset += v.block(block).size();
+          max_size = std::max(max_size, v.block(block).size());
         }
+      owned_elements.set_size(max_size);
 
       Assert (n_elements == added_elements, ExcInternalError());
       Epetra_Map new_map (v.size(), n_elements, &global_ids[0], 0,
@@ -344,6 +380,7 @@ namespace TrilinosWrappers
         {
           v.block(block).trilinos_vector().ExtractCopy (entries, 0);
           entries += v.block(block).local_size();
+          owned_elements.add_indices(v.block(block).locally_owned_elements());
         }
 
       if (import_data == true)
@@ -360,7 +397,6 @@ namespace TrilinosWrappers
 
           last_action = Insert;
         }
-
     }
 
 
@@ -370,11 +406,17 @@ namespace TrilinosWrappers
                         const bool      vector_writable)
     {
       nonlocal_vector.reset();
+      owned_elements = locally_owned_entries;
       if (vector_writable == false)
         {
-          IndexSet parallel_partitioning = locally_owned_entries;
-          parallel_partitioning.add_indices(ghost_entries);
-          reinit(parallel_partitioning, communicator);
+          IndexSet parallel_partitioner = locally_owned_entries;
+          parallel_partitioner.add_indices(ghost_entries);
+          Epetra_Map map = parallel_partitioner.make_trilinos_map (communicator,
+                                                                   true);
+          vector.reset (new Epetra_FEVector(map));
+          has_ghosts = vector->Map().UniqueGIDs()==false;
+
+          last_action = Zero;
         }
       else
         {
@@ -438,6 +480,7 @@ namespace TrilinosWrappers
           vector.reset (new Epetra_FEVector(*v.vector));
           last_action = Zero;
           has_ghosts = v.has_ghosts;
+          owned_elements = v.locally_owned_elements();
         }
 
       if (v.nonlocal_vector.get() != 0)
@@ -519,20 +562,14 @@ namespace TrilinosWrappers
 
   Vector::Vector (const size_type n)
   {
-    last_action = Zero;
-    Epetra_LocalMap map ((TrilinosWrappers::types::int_type)n, 0, Utilities::Trilinos::comm_self());
-    vector.reset (new Epetra_FEVector (map));
+    reinit(n);
   }
 
 
 
   Vector::Vector (const Epetra_Map &input_map)
   {
-    last_action = Zero;
-    Epetra_LocalMap map (n_global_elements(input_map),
-                         input_map.IndexBase(),
-                         input_map.Comm());
-    vector.reset (new Epetra_FEVector(map));
+    reinit(input_map);
   }
 
 
@@ -540,16 +577,7 @@ namespace TrilinosWrappers
   Vector::Vector (const IndexSet &partitioning,
                   const MPI_Comm &communicator)
   {
-    last_action = Zero;
-    Epetra_LocalMap map (static_cast<TrilinosWrappers::types::int_type>(partitioning.size()),
-                         0,
-#ifdef DEAL_II_WITH_MPI
-                         Epetra_MpiComm(communicator));
-#else
-                         Epetra_SerialComm());
-    (void)communicator;
-#endif
-    vector.reset (new Epetra_FEVector(map));
+    reinit (partitioning, communicator);
   }
 
 
@@ -569,7 +597,6 @@ namespace TrilinosWrappers
       }
     else
       reinit (v, false, true);
-
   }
 
 
@@ -578,6 +605,9 @@ namespace TrilinosWrappers
   Vector::reinit (const size_type n,
                   const bool    /*omit_zeroing_entries*/)
   {
+    owned_elements.clear();
+    owned_elements.set_size(n);
+    owned_elements.add_range(0, n);
     Epetra_LocalMap map ((TrilinosWrappers::types::int_type)n, 0,
                          Utilities::Trilinos::comm_self());
     vector.reset (new Epetra_FEVector (map));
@@ -594,7 +624,9 @@ namespace TrilinosWrappers
     Epetra_LocalMap map (n_global_elements(input_map),
                          input_map.IndexBase(),
                          input_map.Comm());
-    vector.reset (new Epetra_FEVector (map));
+    vector.reset (new Epetra_FEVector(map));
+    owned_elements.set_size(n_global_elements(input_map));
+    owned_elements.add_range(0, n_global_elements(input_map));
 
     last_action = Zero;
   }
@@ -617,6 +649,7 @@ namespace TrilinosWrappers
     vector.reset (new Epetra_FEVector(map));
 
     last_action = Zero;
+    owned_elements = partitioning;
   }
 
 
@@ -650,6 +683,7 @@ namespace TrilinosWrappers
                                  v.vector->Map().IndexBase(),
                                  v.vector->Comm());
             vector.reset (new Epetra_FEVector(map));
+            owned_elements = v.locally_owned_elements();
           }
         else if (omit_zeroing_entries)
           {
@@ -723,6 +757,7 @@ namespace TrilinosWrappers
                              v.vector->Map().IndexBase(),
                              v.vector->Comm());
         vector.reset (new Epetra_FEVector(map));
+        owned_elements = v.locally_owned_elements();
       }
 
     const int ierr = vector->Update(1.0, *v.vector, 0.0);
diff --git a/source/lac/trilinos_vector_base.cc b/source/lac/trilinos_vector_base.cc
index fd05084a02..1d7ebe5899 100644
--- a/source/lac/trilinos_vector_base.cc
+++ b/source/lac/trilinos_vector_base.cc
@@ -101,7 +101,10 @@ namespace TrilinosWrappers
     compressed (true),
     has_ghosts  (v.has_ghosts),
     vector(new Epetra_FEVector(*v.vector))
-  {}
+  {
+    owned_elements.set_size(v.size());
+    owned_elements.add_range(0,v.size());
+  }
 
 
 
@@ -139,6 +142,7 @@ namespace TrilinosWrappers
         last_action = Zero;
         vector.reset (new Epetra_FEVector(*v.vector));
         has_ghosts = v.has_ghosts;
+        owned_elements = v.locally_owned_elements();
       }
     else
       {
diff --git a/tests/mpi/trilinos_locally_owned_elements_02.cc b/tests/mpi/trilinos_locally_owned_elements_02.cc
new file mode 100644
index 0000000000..7cd1ace7dc
--- /dev/null
+++ b/tests/mpi/trilinos_locally_owned_elements_02.cc
@@ -0,0 +1,88 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 2016 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// 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>
+
+
+void test ()
+{
+  const int n_proc = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  const int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  //All processes should own 10 entries
+  const int entries_per_process = 10;
+
+  IndexSet locally_owned(entries_per_process*n_proc);
+  const int begin_index = my_id*entries_per_process;
+  const int end_index = (my_id+1)*entries_per_process;
+  locally_owned.add_range(begin_index, end_index);
+
+  IndexSet locally_relevant(entries_per_process*n_proc);
+  const int local_begin = std::max(0, begin_index-entries_per_process/2);
+  const int local_end = entries_per_process*n_proc;
+  locally_relevant.add_range (local_begin, local_end);
+
+  TrilinosWrappers::MPI::Vector ghosted, distributed;
+  distributed.reinit(locally_owned, MPI_COMM_WORLD);
+  ghosted.reinit (locally_owned, locally_relevant, MPI_COMM_WORLD);
+
+  IndexSet locally_owned_elements_distributed = distributed.locally_owned_elements();
+  IndexSet locally_owned_elements_ghosted = ghosted.locally_owned_elements();
+
+  const types::global_dof_index local_range_begin_ghosted = ghosted.local_range().first;
+  const types::global_dof_index local_range_end_ghosted = ghosted.local_range().second;
+
+  const types::global_dof_index local_range_begin_distributed = distributed.local_range().first;
+  const types::global_dof_index local_range_end_distributed = distributed.local_range().second;
+
+  deallog << "locally_owned_elements_distributed: ";
+  locally_owned_elements_distributed.print(deallog);
+  deallog << "locally_owned_elements_ghosted: ";
+  locally_owned_elements_ghosted.print(deallog);
+  deallog << "local_range_begin_ghosted: "
+          << local_range_begin_ghosted << std::endl;
+  deallog << "local_range_end_ghosted: "
+          << local_range_end_ghosted << std::endl;
+  deallog << "local_range_begin_distributed: "
+          << local_range_begin_distributed << std::endl;
+  deallog << "local_range_end_distributed: "
+          << local_range_end_distributed << std::endl;
+
+  AssertThrow (locally_owned_elements_distributed == locally_owned_elements_ghosted,
+               ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+
+  MPILogInitAll log;
+
+  test();
+}
diff --git a/tests/mpi/trilinos_locally_owned_elements_02.with_trilinos=true.mpirun=2.output b/tests/mpi/trilinos_locally_owned_elements_02.with_trilinos=true.mpirun=2.output
new file mode 100644
index 0000000000..ce5508c5b9
--- /dev/null
+++ b/tests/mpi/trilinos_locally_owned_elements_02.with_trilinos=true.mpirun=2.output
@@ -0,0 +1,17 @@
+
+DEAL:0::locally_owned_elements_distributed: {[0,9]}
+DEAL:0::locally_owned_elements_ghosted: {[0,9]}
+DEAL:0::local_range_begin_ghosted: 0
+DEAL:0::local_range_end_ghosted: 20
+DEAL:0::local_range_begin_distributed: 0
+DEAL:0::local_range_end_distributed: 10
+DEAL:0::OK
+
+DEAL:1::locally_owned_elements_distributed: {[10,19]}
+DEAL:1::locally_owned_elements_ghosted: {[10,19]}
+DEAL:1::local_range_begin_ghosted: 5
+DEAL:1::local_range_end_ghosted: 20
+DEAL:1::local_range_begin_distributed: 10
+DEAL:1::local_range_end_distributed: 20
+DEAL:1::OK
+