SF_network.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#ifndef _SF_NETWORK_H
#define _SF_NETWORK_H

#include <mpi.h>

#include "SF_container.h"
#include "SF_globals.h"
#include "SF_vector.h"

namespace SF {

template<class T, class S>
inline void MPI_Exchange(commgraph<T> & grph,
                         vector<S> & send,
                         vector<S> & recv,
                         MPI_Comm comm)
{
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);

  size_t numsend = 0, numrecv = 0;
  for(int i=0; i<size; i++) {
    if(grph.scnt[i]) numsend++;
    if(grph.rcnt[i]) numrecv++;
  }

  vector<MPI_Request> sreq(numsend), rreq(numrecv);
  vector<MPI_Status> stat(numsend > numrecv ? numsend : numrecv);

  size_t sidx = 0, ridx = 0;
  for(int i=0; i<size; i++) {
    if(grph.rcnt[i])
      MPI_Irecv(recv.data() + grph.rdsp[i], grph.rcnt[i]*sizeof(S), MPI_BYTE, i, SF_MPITAG, comm, rreq.data()+ridx++);
    if(grph.scnt[i])
      MPI_Isend(send.data() + grph.sdsp[i], grph.scnt[i]*sizeof(S), MPI_BYTE, i, SF_MPITAG, comm, sreq.data()+sidx++);
  }

  if( (ridx != numrecv) || (sidx != numsend) )
    fprintf(stderr, "Rank %d: MPI_Exchange error! \n", rank);

  MPI_Waitall(ridx, rreq.data(), stat.data());
  MPI_Waitall(sidx, sreq.data(), stat.data());
}


template<class T>
inline T parallel_fallback_value(const vector<T> & vec, MPI_Comm comm)
{
  int rank, size;
  MPI_Comm_rank(comm, &rank); MPI_Comm_size(comm, &size);

  vector<size_t> vec_sizes(size);
  size_t lsize = vec.size();

  MPI_Allgather(&lsize, sizeof(size_t), MPI_BYTE, vec_sizes.data(), sizeof(size_t), MPI_BYTE, comm);

  int fallback_rank = 0;
  while(fallback_rank < size && vec_sizes[fallback_rank] == 0) fallback_rank++;

  assert(vec_sizes[fallback_rank]);

  T fallback_value = T();
  if(rank == fallback_rank)
    fallback_value = vec[0];

  MPI_Bcast(&fallback_value, sizeof(T), MPI_BYTE, fallback_rank, comm);

  return fallback_value;
}


template<class T>
inline T global_min(const vector<T> & vec, MPI_Comm comm)
{
  int size;
  MPI_Comm_size(comm, &size);

  int zero_size_err = 0;
  if(vec.size() == 0) zero_size_err++;
  MPI_Allreduce(MPI_IN_PLACE, &zero_size_err, 1, MPI_INT, MPI_SUM, comm);

  T fb = T();
  if(zero_size_err)
    fb = parallel_fallback_value(vec, comm);

  vector<T> mins(size);

  T lmin = vec.size() ? *std::min_element(vec.begin(), vec.end()) : fb;
  MPI_Allgather(&lmin, sizeof(T), MPI_BYTE, mins.data(), sizeof(T), MPI_BYTE, comm);
  T gmin = *std::min_element(mins.begin(), mins.end());

  return gmin;
}
template<class T>
inline T global_max(const vector<T> & vec, MPI_Comm comm)
{
  int size;
  MPI_Comm_size(comm, &size);

  int zero_size_err = 0;
  if(vec.size() == 0) zero_size_err++;
  MPI_Allreduce(MPI_IN_PLACE, &zero_size_err, 1, MPI_INT, MPI_SUM, comm);

  T fb = T();
  if(zero_size_err)
    fb = parallel_fallback_value(vec, comm);

  vector<T> max(size);

  T lmax = vec.size() ? *std::max_element(vec.begin(), vec.end()) : fb;
  MPI_Allgather(&lmax, sizeof(T), MPI_BYTE, max.data(), sizeof(T), MPI_BYTE, comm);
  T gmax = *std::max_element(max.begin(), max.end());

  return gmax;
}
template<class T>
inline T global_max(const T val, MPI_Comm comm)
{
  int size;
  MPI_Comm_size(comm, &size);
  vector<T> max(size);

  MPI_Allgather(&val, sizeof(T), MPI_BYTE, max.data(), sizeof(T), MPI_BYTE, comm);
  T gmax = *std::max_element(max.begin(), max.end());

  return gmax;
}

template<class T> inline
void layout_from_count(const T count, vector<T> & layout, MPI_Comm comm)
{
  int size;
  MPI_Comm_size(comm, &size);
  vector<T> cnt(size);

  MPI_Allgather(&count, sizeof(T), MPI_BYTE, cnt.data(), sizeof(T), MPI_BYTE, comm);
  dsp_from_cnt(cnt, layout);
}

template<class T> inline
T index_in_layout(const T idx, const vector<T> & layout)
{
  T li = 0, end = layout.size() - 1;

  while(li < end && !(layout[li] <= idx && layout[li+1] > idx))
    li++;

  return li;
}


template<class T> inline
void make_global(const vector<T> & vec, vector<T> & out, MPI_Comm comm)
{
  int vecsize = vec.size()*sizeof(T);
  vector<int> cnt, dsp;
  layout_from_count(vecsize, dsp, comm);
  cnt_from_dsp(dsp, cnt);

  out.resize(sum(cnt)/sizeof(T));
  MPI_Allgatherv(vec.data(), vecsize, MPI_BYTE, out.data(), cnt.data(), dsp.data(), MPI_BYTE, comm);
}

template<class T> inline
void make_global(vector<T> & vec, MPI_Comm comm)
{
  vector<T> out;
  make_global(vec, out, comm);
  vec = out;
}

}
#endif