Mercurial > fem-fenics-eugenio
changeset 256:8fe68d94ab76
Add parallel assembly of matrices and vectors
author | Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> |
---|---|
date | Wed, 30 Jul 2014 20:59:37 +0200 |
parents | 072aee55bb75 |
children | fb67b636616f |
files | src/Function.cc src/Makefile.in src/PETSc_factory.cc src/PETSc_factory.h |
diffstat | 4 files changed, 257 insertions(+), 42 deletions(-) [+] |
line wrap: on
line diff
--- 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; +}
--- 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)'
--- 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; + }
--- 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