librheolef-dev 6.7-6 / usr / include / rheolef / mpi_assembly_begin.h

#ifndef _RHEO_MPI_ASSEMBLY_BEGIN_H
#define _RHEO_MPI_ASSEMBLY_BEGIN_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================

namespace rheolef {

/*F:
NAME: mpi_assembly_begin -- for array or matrix (@PACKAGE@ @VERSION@)
DESCRIPTION:
  Start a dense array or a sparse matrix assembly.
COMPLEXITY: 
  Assume that stash has indexes in increasing order.
  cpu complexity : O(stash.size + nproc)
  memory complexity : O(stash.size + nproc)
  msg size complexity : O(stash.size + nproc)
  
  When assembling a finite element matrix, the stash.size is at boundaries
  of the mesh partition, and stash.size = O(nrow^alpha), where alpha=(d-1)/d
  and d=2,3. Using nproc <= O(nrow^alpha) is performant.
NOTE:
  The stash may be sorted by increasing nows and column.

  Inspirated from Petsc (petsc/src/vec/vec/impls/mpi/pdvec.c), here with
  a pre-sorted stash, thanks to the stl::map data structure.
AUTHORS:
    LMC-IMAG, 38041 Grenoble cedex 9, France
    | Pierre.Saramito@imag.fr
DATE:   23 march 1999
END:
*/

//<mpi_assembly_begin:
template <
    class Stash,
    class Message,
    class InputIterator>
typename Stash::size_type
mpi_assembly_begin (
 // input:
    const Stash&                       stash,
    InputIterator                      first_stash_idx,	// wrapper in shash
    InputIterator                      last_stash_idx,	
    const distributor&                 ownership,
 // ouput:
    Message&                           receive,		// buffer
    Message&                           send)		// buffer
{
    typedef typename Stash::size_type      size_type;

    mpi::communicator     comm    = ownership.comm();
    size_type             my_proc = ownership.process();
    size_type             nproc   = ownership.n_process();
    distributor::tag_type tag     = distributor::get_new_tag();

    // ----------------------------------------------------------------
    // 1) count the messages contained in stash by process id
    // ----------------------------------------------------------------
    // assume that stash elements are sorted by increasing stash_idx (e.g. stash = stl::map)
    // cpu complexity = O(stash.size + nproc)
    // mem complexity = O(stash.size + nproc)
    std::vector<size_type> msg_size(nproc, 0);
    std::vector<size_type> msg_mark(nproc, 0);
    size_type send_nproc = 0;
    {
      size_type iproc = 0;
      size_type i = 0;
      for (InputIterator iter_idx = first_stash_idx; iter_idx != last_stash_idx; iter_idx++, i++) {
        for (; iproc < nproc; iproc++) {
          if (*iter_idx >= ownership[iproc] && *iter_idx < ownership[iproc+1]) {
            msg_size [iproc]++;
            if (!msg_mark[iproc]) {
               msg_mark[iproc] = 1;
               send_nproc++;
            }
            break;
          }
        }
        assert_macro (iproc != nproc, "bad stash data: index "<<*iter_idx<<" out of range [0:"<<ownership[nproc]<<"[");
      }
    } // end block
    // ----------------------------------------------------------------
    // 2) avoid to send message to my-proc in counting
    // ----------------------------------------------------------------
    if (msg_size [my_proc] != 0) {
	msg_size [my_proc] = 0;
	msg_mark [my_proc] = 0;
	send_nproc--;
    }
    // ----------------------------------------------------------------
    // 3) compute number of messages to be send to my_proc
    // ----------------------------------------------------------------
    // msg complexity : O(nproc) or O(log(nproc)), depending on reduce implementation 
    std::vector<size_type> work (nproc);
    mpi::all_reduce (
	comm, 
        msg_mark.begin().operator->(),
	nproc,
	work.begin().operator->(),
	std::plus<size_type>());
    size_type receive_nproc = work [my_proc];
    // ----------------------------------------------------------------
    // 4) compute messages max size to be send to my_proc
    // ----------------------------------------------------------------
    // msg complexity : O(nproc) or O(log(nproc)), depending on reduce implementation 
    mpi::all_reduce (
        comm,
        msg_size.begin().operator->(),
        nproc,
	work.begin().operator->(),
        mpi::maximum<size_type>());
    size_type receive_max_size = work [my_proc];
    // ----------------------------------------------------------------
    // 5) post receive: exchange the buffer adresses between processes
    // ----------------------------------------------------------------
    // Each message will consist of ordered pairs (global index,value).
    // since we don't know how long each indiidual message is,
    // we allocate the largest : receive_nproc*receive_max_size
    // potentially, this is a lot of wasted space
    // TODO: how to optimize the receive.data buffer ?
    // cpu complexity : O(nproc)
    // mem complexity : O(nproc*(stash.size/nproc)) = O(stash.size), worst case ?
    // msg complexity : O(nproc) 
    receive.data.resize (receive_nproc*receive_max_size);
    for (size_t i_receive = 0; i_receive < receive_nproc; i_receive++) {
      mpi::request i_req = comm.irecv (
	  mpi::any_source,
	  tag,
          receive.data.begin().operator->() + i_receive*receive_max_size,
	  receive_max_size);
      receive.waits.push_back (std::make_pair(i_receive, i_req));
    }
    // ----------------------------------------------------------------
    // 6) copy stash in send buffer
    // ----------------------------------------------------------------
    // since the stash is sorted by increasing order => simple copy
    // cpu complexity : O(stash.size)
    // mem complexity : O(stash.size)
    send.data.resize (stash.size());
    copy (stash.begin(), stash.end(), send.data.begin());
    // ---------------------------------------------------------------------------
    // 7) do send
    // ---------------------------------------------------------------------------
    // cpu complexity : O(nproc)
    // mem complexity : O(send_nproc) \approx O(nproc), worst case
    // msg complexity : O(stash.size)
    send.waits.resize(send_nproc);
    {
      size_type i_send = 0;
      size_type i_start = 0;
      for (size_type iproc = 0; iproc < nproc; iproc++) {
        size_type i_msg_size = msg_size[iproc];
        if (i_msg_size == 0) continue;
        mpi::request i_req = comm.isend (
	    iproc, 
	    tag, 
	    send.data.begin().operator->() + i_start, 
	    i_msg_size);
        send.waits.push_back(std::make_pair(i_send,i_req));
        i_send++;
        i_start += i_msg_size;
      }
    } // end block
    return receive_max_size;
}
//>mpi_assembly_begin:
} // namespace rheolef
#endif //_RHEO_MPI_ASSEMBLY_BEGIN_H