# HG changeset patch
# User Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
# Date 1406746777 -7200
# Node ID 8fe68d94ab76cc7a520b6e84a726d16bbc621399
# Parent  072aee55bb7575b334f96418c2b6c0721c6adcd8
Add parallel assembly of matrices and vectors

diff -r 072aee55bb75 -r 8fe68d94ab76 src/Function.cc
--- a/src/Function.cc	Wed Jul 30 18:09:13 2014 +0200
+++ b/src/Function.cc	Wed Jul 30 20:59:37 2014 +0200
@@ -19,6 +19,13 @@
 #include "functionspace.h"
 #include "function.h"
 #include "dolfin_compat.h"
+#ifdef LATEST_DOLFIN
+#include <boost/mpi.hpp>
+#include <boost/serialization/vector.hpp>
+#include <boost/serialization/utility.hpp>
+#endif
+
+void copy_vector (Array <double> const&, SHARED_PTR <dolfin::GenericVector>);
 
 DEFUN_DLD (Function, args, , "-*- texinfo -*-\n\
 @deftypefn {Function File} {[@var{func}]} = \
@@ -90,20 +97,11 @@
                        and of the vector must agree");
               else
                 {
-#ifdef LATEST_DOLFIN
-                  dolfin::Vector du (MPI_COMM_WORLD, myu.length ());
-#else
-                  dolfin::Vector du(myu.length ());
-#endif
-                  for (std::size_t i = 0; i < myu.length (); ++i)
-                    du.setitem (i, myu(i));
+                  dolfin::Function aux (V);
+                  copy_vector (myu, aux.vector ());
 
-                  SHARED_PTR <dolfin::Vector> 
-                    uu (new dolfin::Vector(du));
-
-                  SHARED_PTR <const dolfin::Function>
-                    u (new dolfin::Function(V, uu));
-
+                  SHARED_PTR <dolfin::Function const>
+                    u (new dolfin::Function (aux));
                   retval = new function (str, u);
                }
             }
@@ -140,3 +138,64 @@
     }
   return retval;
 }
+
+void
+copy_vector (Array <double> const& arr, SHARED_PTR <dolfin::GenericVector> vec)
+{
+  unsigned const size =
+#ifdef LATEST_DOLFIN
+    dolfin::MPI::size (MPI_COMM_WORLD);
+#else
+    dolfin::MPI::num_processes ();
+#endif
+
+  if (size > 1)
+    {
+      std::vector <double> locvals;
+      boost::mpi::communicator world;
+      unsigned const rank = world.rank ();
+
+      std::pair <std::size_t, std::size_t> const loc_range =
+        vec->local_range ();
+
+      for (unsigned p = 1; p < size; ++p)
+        if (rank == p)
+          world.send (0, p, loc_range);
+
+      if (rank == 0)
+        {
+          for (unsigned p = 0; p < size; ++p)
+            {
+              std::pair <std::size_t, std::size_t> range;
+              if (p == 0)
+                range = loc_range;
+              else
+                world.recv (p, p, range);
+
+              std::vector <double> entries;
+              for (std::size_t i = range.first; i < range.second; ++i)
+                entries.push_back (arr(i));
+
+              if (p == 0)
+                locvals = entries;
+              else
+                world.send (p, p, entries);
+            }
+        }
+
+      for (unsigned p = 1; p < size; ++p)
+        {
+          if (rank == p)
+            world.recv (0, p, locvals);
+        }
+
+      vec->set_local (locvals);
+    }
+  else
+    for (std::size_t i = 0; i < arr.length (); ++i)
+      vec->setitem (i, arr(i));
+
+  vec->apply ("insert");
+
+  return;
+}
diff -r 072aee55bb75 -r 8fe68d94ab76 src/Makefile.in
--- a/src/Makefile.in	Wed Jul 30 18:09:13 2014 +0200
+++ b/src/Makefile.in	Wed Jul 30 20:59:37 2014 +0200
@@ -24,6 +24,7 @@
 LIBS=$(patsubst %, "%", $(LIBS_RAW))
 ifeq (1.4.0, $(findstring 1.4.0, $(DOLFIN_CPPFLAGS)))
   CPPFLAGS+=-DLATEST_DOLFIN
+  LIBS+=-lboost_mpi -lboost_serialization
 endif
 CPPFLAGS+=$(DOLFIN_CPPFLAGS)
 CPPFLAGS:='$(CPPFLAGS)'
diff -r 072aee55bb75 -r 8fe68d94ab76 src/PETSc_factory.cc
--- a/src/PETSc_factory.cc	Wed Jul 30 18:09:13 2014 +0200
+++ b/src/PETSc_factory.cc	Wed Jul 30 20:59:37 2014 +0200
@@ -16,6 +16,10 @@
 */
 
 #include "PETSc_factory.h"
+#ifdef LATEST_DOLFIN
+#include <boost/mpi.hpp>
+#include <boost/serialization/vector.hpp>
+#endif
 
 femfenics_base_factory const&
 PETSc_factory::instance (void)
@@ -29,25 +33,17 @@
   {
     octave_value retval;
 
-    std::size_t nr = A.size (0), nc = A.size (1);
-
-    dim_vector dims (0, 1);
-    Array <octave_idx_type> 
-      ridx (dims, 0), 
-      cidx (dims, 0);
-    Array<double> data (dims, 0);
+    unsigned const size =
+#ifdef LATEST_DOLFIN
+      dolfin::MPI::size (MPI_COMM_WORLD);
+#else
+      dolfin::MPI::num_processes ();
+#endif
 
-    for (std::size_t i = 0; i < nr; ++i)
-      {
-        std::vector <double> data_tmp;
-        std::vector <std::size_t> cidx_tmp;
-        A.getrow (i, cidx_tmp, data_tmp);
-        std::vector <std::size_t> ridx_tmp (cidx_tmp.size (), i);
-        add_to_arrays (ridx, cidx, data, ridx_tmp, cidx_tmp, data_tmp);
-      }
-
-    SparseMatrix sm (data, ridx, cidx, nr, nc);
-    retval = sm;
+    if (size > 1)
+      retval = do_matrix_parallel (A);
+    else
+      retval = do_matrix_serial (A);
 
     return retval;
   }
@@ -57,16 +53,17 @@
   {
     octave_value retval;
 
-    dim_vector dims;
-    dims.resize (2);
-    dims(0) = b.size ();
-    dims(1) = 1;
-    Array <double> myb (dims);
+    unsigned const size =
+#ifdef LATEST_DOLFIN
+      dolfin::MPI::size (MPI_COMM_WORLD);
+#else
+      dolfin::MPI::num_processes ();
+#endif
 
-    for (std::size_t i = 0; i < b.size (); ++i)
-      myb.xelem (i) = b[i];
-
-    retval = myb;
+    if (size > 1)
+      retval = do_vector_parallel (b);
+    else
+      retval = do_vector_serial (b);
 
     return retval;
   }
@@ -107,4 +104,157 @@
       }
 
     return;
-  }
\ No newline at end of file
+  }
+
+octave_value
+PETSc_factory::do_matrix_serial (dolfin::Matrix const& A) const
+  {
+    octave_value retval;
+
+    std::size_t nr = A.size (0), nc = A.size (1);
+
+    dim_vector dims (0, 1);
+    Array <octave_idx_type> 
+      ridx (dims, 0), 
+      cidx (dims, 0);
+    Array<double> data (dims, 0);
+
+    for (std::size_t i = 0; i < nr; ++i)
+      {
+        std::vector <double> data_tmp;
+        std::vector <std::size_t> cidx_tmp;
+        A.getrow (i, cidx_tmp, data_tmp);
+        std::vector <std::size_t> ridx_tmp (cidx_tmp.size (), i);
+        add_to_arrays (ridx, cidx, data, ridx_tmp, cidx_tmp, data_tmp);
+      }
+
+    SparseMatrix sm (data, ridx, cidx, nr, nc);
+    retval = sm;
+
+    return retval;
+  }
+
+octave_value
+PETSc_factory::do_matrix_parallel (dolfin::Matrix const& A) const
+  {
+    octave_value retval;
+
+    std::vector <std::size_t> rows, cols;
+    std::vector <double> vals;
+    std::size_t const nrows = A.size (0), ncols = A.size (1);
+
+    std::pair <std::size_t, std::size_t> const range = A.local_range (0);
+    for (std::size_t i = range.first; i < range.second; ++i)
+      {
+        std::vector <std::size_t> cols_tmp;
+        std::vector <double> vals_tmp;
+        A.getrow (i, cols_tmp, vals_tmp);
+        cols.insert (cols.end (), cols_tmp.begin (), cols_tmp.end ());
+        vals.insert (vals.end (), vals_tmp.begin (), vals_tmp.end ());
+        for (std::size_t j = 0; j < cols_tmp.size (); ++j)
+          rows.push_back (i);
+      }
+
+    boost::mpi::communicator world;
+
+    unsigned const size = world.size ();
+    unsigned const rank = world.rank ();
+
+    for (unsigned p = 1; p < size; ++p)
+      {
+        if (rank == p)
+          {
+            unsigned const tag = 3*(p-1);
+            world.send (0, tag, rows);
+            world.send (0, tag+1, cols);
+            world.send (0, tag+2, vals);
+          }
+      }
+
+    if (rank == 0)
+      {
+        dim_vector dims (0, 1);
+        Array <octave_idx_type> ridx (dims, 0), cidx (dims, 0);
+        Array <double> data (dims, 0);
+        add_to_arrays (ridx, cidx, data, rows, cols, vals);
+
+        for (unsigned p = 1; p < size; ++p)
+          {
+            std::vector <double> new_vals;
+            std::vector <std::size_t> new_rows, new_cols;
+            unsigned const tag = 3*(p-1);
+
+            world.recv (p, tag, new_rows);
+            world.recv (p, tag+1, new_cols);
+            world.recv (p, tag+2, new_vals);
+
+            add_to_arrays (ridx, cidx, data, new_rows, new_cols, new_vals);
+          }
+
+        SparseMatrix sm (data, ridx, cidx, nrows, ncols);
+        retval = sm;
+      }
+    else
+      // Return an identity matrix just to avoid warnings or errors
+      {
+        dim_vector dims (nrows, 1);
+        Array <octave_idx_type> ridx (dims), cidx (dims);
+        Array <double> data (dims, 1);
+
+        for (std::size_t i = 0; i < nrows; ++i)
+          {
+            ridx.xelem (i) = i;
+            cidx.xelem (i) = i;
+          }
+
+        SparseMatrix sm (data, ridx, cidx, nrows, ncols);
+        retval = sm;
+      }
+
+    return retval;
+  }
+
+octave_value
+PETSc_factory::do_vector_serial (dolfin::Vector const& b) const
+  {
+    octave_value retval;
+
+    dim_vector dims;
+    dims.resize (2);
+    dims(0) = b.size ();
+    dims(1) = 1;
+    Array <double> myb (dims);
+
+    for (std::size_t i = 0; i < b.size (); ++i)
+      myb.xelem (i) = b[i];
+
+    retval = myb;
+
+    return retval;
+  }
+
+octave_value
+PETSc_factory::do_vector_parallel (dolfin::Vector const& b) const
+  {
+    octave_value retval;
+
+    std::size_t const length = b.size ();
+    dim_vector dims (length, 1);
+    Array <double> vec (dims, 0);
+
+    std::vector <double> entries;
+    b.gather_on_zero (entries);
+
+#ifdef LATEST_DOLFIN
+    if (dolfin::MPI::rank (MPI_COMM_WORLD) == 0)
+#else
+    if (dolfin::MPI::process_number () == 0)
+#endif
+      for (std::size_t i = 0; i < length; ++i)
+        vec.xelem (i) = entries[i];
+
+    // On nodes other than zero, return an all zero vector of the right length
+    retval = vec;
+
+    return retval;
+  }
diff -r 072aee55bb75 -r 8fe68d94ab76 src/PETSc_factory.h
--- a/src/PETSc_factory.h	Wed Jul 30 18:09:13 2014 +0200
+++ b/src/PETSc_factory.h	Wed Jul 30 20:59:37 2014 +0200
@@ -43,6 +43,11 @@
                              std::vector <std::size_t> const&,
                              std::vector <std::size_t> const&,
                              std::vector <double> const&);
+
+  octave_value do_matrix_serial (dolfin::Matrix const&) const;
+  octave_value do_matrix_parallel (dolfin::Matrix const&) const;
+  octave_value do_vector_serial (dolfin::Vector const&) const;
+  octave_value do_vector_parallel (dolfin::Vector const&) const;
 };
 
-#endif
\ No newline at end of file
+#endif