bgpdma.cpp

Go to the documentation of this file.

// included by netpar.cpp
 
/*
Overall exchange strategy
 
When a cell spikes, it immediately does a DCMF_Multicast of
(int gid, double spiketime) to all the target machines that have
cells that need to receive this spike by spiketime + delay
I'd like to cycle through a list of mconfig.nconnections so that
I don't have to wait for my single connection to complete the previous
broadcast when there is a high density of generated spikes but I need
to take care of my bus error issues first.
 
In order to minimize the number of nrnmpi_bgp_conserve tests
(and potentially abandon them altogether if I can ever guarantee
that exchange time is less than half the computation time), I divide the
minimum delay integration intervals into two equal subintervals.
So if a spike is generated in an even subinterval, I do not have
to include it in the conservation check until the end of the next even
subinterval.
 
When a spike is received (DMA interrupt) it is placed in even or odd
buffers (depending on whether the coded gid is positive or negative)
 
At the end of a computation subinterval the even or odd buffer spikes
are enqueued in the priority queue after checking that the number
of spikes sent is equal to the number of spikes sent.
*/
 
extern "C" {
 
extern IvocVect* vector_arg(int);
extern void vector_resize(IvocVect*, int);
 
}  // extern "C"
extern void (*nrntimeout_call)();
 
 
// The initial idea behind use_phase2_ is to avoid the large overhead of
// initiating a send of the up to 10k list of target hosts when a cell fires.
// I.e. when there are a small number of cells on a processor, this causes
// load balance problems.
// Load balance shuld be better if the send is distributed to a much smaller
// set of targets, which, when they receive the spike, pass it on to a neighbor
// set.
// We expect that TWOPHASE will work best in combination with ENQUEUE=2
// which has the greatest amount of overlap between computation
// and communication.
// Note: the old implementation assumed that input PreSyn did not need
// a BGP_DMASend pointer so used the PreSyn.bgp.srchost_ element to
// help figure out the target_hosts_ list for the output PreSyn.bgp.dma_send_.
// Since bgp.srchost_ is used only for setup, it can be overwritten at phase2
// setup time with a bgp.dma_send_ so as to pass on the spike to the
// phase2 list of target hosts.
 
#if !defined(DCMFTICK)
#define DCMFTICK     0
#define DCMFTIMEBASE 0
#endif
 
static unsigned long long dmasend_time_;
static int n_xtra_cons_check_;
#define MAXNCONS 10
#if MAXNCONS
static int xtra_cons_hist_[MAXNCONS + 1];
#endif
 
// asm/msr.h no longer compiles on my machine.
// only for basic testing of logic when not on blue gene/p
#define USE_RDTSCL 0
 
// only use if careful not to overrun the buffer during a simulation
#if 0 && USE_RDTSCL
#define TBUFSIZE (1 << 15)
#else
#define TBUFSIZE 0
#endif
 
#if TBUFSIZE
static unsigned long tbuf_[TBUFSIZE];
static int itbuf_;
#if USE_RDTSCL  // but can have many accuracy problems on recent cpus.
/* for rdtscl() */
//#include <asm/msr.h>
#define rdtscl(a) a++
static unsigned long t__;
#define TBUF                   \
    {                          \
        rdtscl(t__);           \
        tbuf_[itbuf_++] = t__; \
    }
#else
#define TBUF tbuf_[itbuf_++] = (unsigned long) DCMF_Timebase();
#endif  // not USE_RDTSCLL
#else
#define TBUF /**/
#endif
 
// ENQUEUE 0 means to  BGP_ReceiveBuffer buffer -> PreSyn.send
// ENQUEUE 1 means to BGP_ReceiveBuffer buffer -> psbuf -> PreSyn.send
// ENQUEUE 2 means to BGP_ReceiveBuffer.incoming -> PrySyn.send
// Note that ENQUEUE 2 give more overlap between computation and exchange
// since the enqueuing takes place during computation except for those
// remaining during conservation.
#define ENQUEUE 2
 
#if ENQUEUE == 2
static unsigned long enq2_find_time_;
static unsigned long enq2_enqueue_time_;  // includes enq_find_time_
#endif
 
#define PHASE2BUFFER_SIZE 2048  // power of 2
#define PHASE2BUFFER_MASK (PHASE2BUFFER_SIZE - 1)
struct Phase2Buffer {
    PreSyn* ps;
    double spiketime;
};
 
#include <structpool.h>
 
using SpkPool = Pool<NRNMPI_Spike>;
 
#define BGP_RECEIVEBUFFER_SIZE 10000
class BGP_ReceiveBuffer {
  public:
    BGP_ReceiveBuffer();
    virtual ~BGP_ReceiveBuffer();
    void init(int index);
    void incoming(int gid, double spiketime);
    void enqueue();
    int index_;
    int size_;
    int count_;
    int maxcount_;
    int busy_;
    int nsend_, nrecv_;  // for checking conservation
    int nsend_cell_;     // cells that spiked this interval.
    unsigned long long timebase_;
    NRNMPI_Spike** buffer_;
    SpkPool* pool_;
 
#if ENQUEUE == 1
    void enqueue1();
    void enqueue2();
#endif
    PreSyn** psbuf_;
    void phase2send();
    int phase2_head_;
    int phase2_tail_;
    int phase2_nsend_cell_, phase2_nsend_;
    Phase2Buffer* phase2_buffer_;
};
 
static int use_phase2_;
#define NTARGET_HOSTS_PHASE1 ntarget_hosts_phase1_
 
class BGP_DMASend {
  public:
    BGP_DMASend();
    virtual ~BGP_DMASend();
    void send(int gid, double t);
    int ntarget_hosts_;
    int* target_hosts_;
    NRNMPI_Spike spk_;
    int ntarget_hosts_phase1_;
};
 
class BGP_DMASend_Phase2 {
  public:
    BGP_DMASend_Phase2();
    virtual ~BGP_DMASend_Phase2();
    NRNMPI_Spike spk_;
 
    void send_phase2(int gid, double t, BGP_ReceiveBuffer*);
    int ntarget_hosts_phase2_;
    int* target_hosts_phase2_;
};
 
static BGP_ReceiveBuffer* bgp_receive_buffer[BGP_INTERVAL];
static int current_rbuf, next_rbuf;
#if BGP_INTERVAL == 2
// note that if a spike is supposed to be received by bgp_receive_buffer[1]
// then during transmission its gid is complemented.
#endif
 
BGP_ReceiveBuffer::BGP_ReceiveBuffer() {
    busy_ = 0;
    count_ = 0;
    size_ = BGP_RECEIVEBUFFER_SIZE;
    buffer_ = new NRNMPI_Spike*[size_];
    pool_ = new SpkPool(BGP_RECEIVEBUFFER_SIZE);
    psbuf_ = 0;
#if ENQUEUE == 1
    psbuf_ = new PreSyn*[size_];
#endif
    phase2_buffer_ = new Phase2Buffer[PHASE2BUFFER_SIZE];
    phase2_head_ = phase2_tail_ = 0;
}
BGP_ReceiveBuffer::~BGP_ReceiveBuffer() {
    assert(busy_ == 0);
    for (int i = 0; i < count_; ++i) {
        pool_->hpfree(buffer_[i]);
    }
    delete[] buffer_;
    delete pool_;
    if (psbuf_)
        delete[] psbuf_;
    delete[] phase2_buffer_;
}
void BGP_ReceiveBuffer::init(int index) {
    index_ = index;
    timebase_ = 0;
    nsend_cell_ = nsend_ = nrecv_ = busy_ = maxcount_ = 0;
    for (int i = 0; i < count_; ++i) {
        pool_->hpfree(buffer_[i]);
    }
    count_ = 0;
    phase2_head_ = phase2_tail_ = 0;
    phase2_nsend_cell_ = phase2_nsend_ = 0;
}
void BGP_ReceiveBuffer::incoming(int gid, double spiketime) {
    // printf("%d %p.incoming %g %g %d\n", nrnmpi_myid, this, t, spk->spiketime, spk->gid);
    assert(busy_ == 0);
    busy_ = 1;
    if (count_ >= size_) {
        size_ *= 2;
        NRNMPI_Spike** newbuf = new NRNMPI_Spike*[size_];
        for (int i = 0; i < count_; ++i) {
            newbuf[i] = buffer_[i];
        }
        delete[] buffer_;
        buffer_ = newbuf;
        if (psbuf_) {
            delete[] psbuf_;
            psbuf_ = new PreSyn*[size_];
        }
    }
    NRNMPI_Spike* spk = pool_->alloc();
    spk->gid = gid;
    spk->spiketime = spiketime;
    buffer_[count_++] = spk;
    if (maxcount_ < count_) {
        maxcount_ = count_;
    }
    ++nrecv_;
    busy_ = 0;
}
void BGP_ReceiveBuffer::enqueue() {
    // printf("%d %p.enqueue count=%d t=%g nrecv=%d nsend=%d\n", nrnmpi_myid, this, t, count_,
    // nrecv_, nsend_);
    assert(busy_ == 0);
    busy_ = 1;
#if 1
    for (int i = 0; i < count_; ++i) {
        NRNMPI_Spike* spk = buffer_[i];
#if ENQUEUE == 2
        unsigned long long tb = DCMFTIMEBASE;
#endif
 
        auto iter = gid2in_.find(spk->gid);
        nrn_assert(iter != gid2in_.end());
        PreSyn* ps = iter->second;
        if (use_phase2_ && ps->bgp.dma_send_phase2_) {
            // cannot do directly because busy_;
            // ps->bgp.dma_send_phase2_->send_phase2(spk->gid, spk->spiketime, this);
            Phase2Buffer& pb = phase2_buffer_[phase2_head_++];
            phase2_head_ &= PHASE2BUFFER_MASK;
            assert(phase2_head_ != phase2_tail_);
            pb.ps = ps;
            pb.spiketime = spk->spiketime;
        }
#if ENQUEUE == 2
        enq2_find_time_ += (unsigned long) (DCMFTIMEBASE - tb);
#endif
        ps->send(spk->spiketime, net_cvode_instance, nrn_threads);
        pool_->hpfree(spk);
#if ENQUEUE == 2
        enq2_enqueue_time_ += (unsigned long) (DCMFTIMEBASE - tb);
#endif
    }
#endif
    count_ = 0;
#if ENQUEUE != 2
    nrecv_ = 0;
    nsend_ = 0;
    nsend_cell_ = 0;
#endif
    busy_ = 0;
    phase2send();
}
 
#if ENQUEUE == 1
void BGP_ReceiveBuffer::enqueue1() {
    // printf("%d %lx.enqueue count=%d t=%g nrecv=%d nsend=%d\n", nrnmpi_myid, (long)this, t,
    // count_, nrecv_, nsend_);
    assert(busy_ == 0);
    busy_ = 1;
    for (int i = 0; i < count_; ++i) {
        NRNMPI_Spike* spk = buffer_[i];
        auto iter = gid2in_->find(spk->gid);
        nrn_assert(iter != gid2in_.end()));
        PreSyn* ps = iter->second;
        psbuf_[i] = ps;
        if (use_phase2_ && ps->bgp.dma_send_phase2_) {
            // cannot do directly because busy_;
            // ps->bgp.dma_send_phase2_->send_phase2(spk->gid, spk->spiketime, this);
            Phase2Buffer& pb = phase2_buffer_[phase2_head_++];
            phase2_head_ &= PHASE2BUFFER_MASK;
            assert(phase2_head_ != phase2_tail_);
            pb.ps = ps;
            pb.spiketime = spk->spiketime;
        }
    }
    busy_ = 0;
    phase2send();
}
 
void BGP_ReceiveBuffer::enqueue2() {
    // printf("%d %lx.enqueue count=%d t=%g nrecv=%d nsend=%d\n", nrnmpi_myid, (long)this, t,
    // count_, nrecv_, nsend_);
    assert(busy_ == 0);
    busy_ = 1;
    for (int i = 0; i < count_; ++i) {
        NRNMPI_Spike* spk = buffer_[i];
        PreSyn* ps = psbuf_[i];
        ps->send(spk->spiketime, net_cvode_instance, nrn_threads);
        pool_->hpfree(spk);
    }
    count_ = 0;
    nrecv_ = 0;
    nsend_ = 0;
    nsend_cell_ = 0;
    busy_ = 0;
}
#endif  // ENQUEUE == 1
 
void BGP_ReceiveBuffer::phase2send() {
    while (phase2_head_ != phase2_tail_) {
        Phase2Buffer& pb = phase2_buffer_[phase2_tail_++];
        phase2_tail_ &= PHASE2BUFFER_MASK;
        pb.ps->bgp.dma_send_phase2_->send_phase2(pb.ps->gid_, pb.spiketime, this);
    }
}
 
// number of DCMF_Multicast_t to cycle through when not using recordreplay
#define NSEND 10
 
static int max_ntarget_host;
// For one phase sending, max_multisend_targets is max_ntarget_host.
// For two phase sending, it is the maximum of all the
// ntarget_hosts_phase1 and ntarget_hosts_phase2.
static int max_multisend_targets;
 
double nrn_bgp_receive_time(int type) {  // and others
    double rt = 0.;
    switch (type) {
    case 2:  // in msend_recv
        if (!use_bgpdma_) {
            return rt;
        }
        for (int i = 0; i < n_bgp_interval; ++i) {
            rt += bgp_receive_buffer[i]->timebase_ * DCMFTICK;
        }
        break;
    case 3:  // in BGP_DMAsend::send
        if (!use_bgpdma_) {
            return rt;
        }
        rt = dmasend_time_ * DCMFTICK;
        break;
    case 4:  // number of extra conservation checks
        rt = double(n_xtra_cons_check_);
        // and if there is second vector arg then also return the histogram
#if MAXNCONS
        if (ifarg(2) && use_bgpdma_) {
            IvocVect* vec = vector_arg(2);
            vector_resize(vec, MAXNCONS + 1);
            for (int i = 0; i <= MAXNCONS; ++i) {
                vector_vec(vec)[i] = double(xtra_cons_hist_[i]);
            }
        }
#endif  // MAXNCONS
#if TBUFSIZE
        if (ifarg(3)) {
            IvocVect* vec = vector_arg(3);
            vector_resize(vec, itbuf_ + 1);
            for (int i = 0; i <= itbuf_; ++i) {
                vector_vec(vec)[i] = double(tbuf_[i]);
            }
            vector_vec(vec)[itbuf_] = DCMFTICK;
        }
#endif
        break;
    case 8:  // exchange method properties
             // bit 0: 0 allgather, 1 multisend (MPI_ISend)
             // bit 1: unused, legacy
             // bit 2: n_bgp_interval, 0 means one interval, 1 means 2
             // bit 3: number of phases, 0 means 1 phase, 1 means 2
             // bit 4: unused (1 used to mean althash used)
             // bit 5: 1 means enqueue separated into two parts for timeing
    {
        int method = use_bgpdma_ ? 1 : 0;
        int p = method + 4 * (n_bgp_interval == 2 ? 1 : 0) + 8 * use_phase2_ +
                16 * (0)  // no hash selection, just std::unordered_map
                + 32 * ENQUEUE;
        rt = double(p);
    } break;
    case 12:  // greatest length multisend
    {
        rt = double(max_multisend_targets);
        break;
    }
    }
    return rt;
}
 
extern void nrnmpi_bgp_comm();
extern void nrnmpi_bgp_multisend(NRNMPI_Spike*, int, int*);
extern int nrnmpi_bgp_single_advance(NRNMPI_Spike*);
extern int nrnmpi_bgp_conserve(int nsend, int nrecv);
 
static void bgp_dma_init() {
    for (int i = 0; i < n_bgp_interval; ++i) {
        bgp_receive_buffer[i]->init(i);
    }
    current_rbuf = 0;
    next_rbuf = n_bgp_interval - 1;
#if TBUFSIZE
    itbuf_ = 0;
#endif
    dmasend_time_ = 0;
#if ENQUEUE == 2
    enq2_find_time_ = enq2_enqueue_time_ = 0;
#endif
    n_xtra_cons_check_ = 0;
#if MAXNCONS
    for (int i = 0; i <= MAXNCONS; ++i) {
        xtra_cons_hist_[i] = 0;
    }
#endif  // MAXNCONS
}
 
static int bgp_advance() {
    NRNMPI_Spike spk;
    int i = 0;
    while (nrnmpi_bgp_single_advance(&spk)) {
        i += 1;
        int j = 0;
#if BGP_INTERVAL == 2
        if (spk.gid < 0) {
            spk.gid = ~spk.gid;
            j = 1;
        }
#endif
        bgp_receive_buffer[j]->incoming(spk.gid, spk.spiketime);
    }
    nrecv_ += i;
    return i;
}
 
#if BGPDMA
void nrnbgp_messager_advance() {
#if BGPDMA & 1
    if (use_bgpdma_) {
        bgp_advance();
    }
#endif
#if ENQUEUE == 2
    bgp_receive_buffer[current_rbuf]->enqueue();
#endif
}
#endif
 
BGP_DMASend::BGP_DMASend() {
    ntarget_hosts_ = 0;
    target_hosts_ = NULL;
    ntarget_hosts_phase1_ = 0;
}
 
BGP_DMASend::~BGP_DMASend() {
    if (target_hosts_) {
        delete[] target_hosts_;
    }
}
 
BGP_DMASend_Phase2::BGP_DMASend_Phase2() {
    ntarget_hosts_phase2_ = 0;
    target_hosts_phase2_ = 0;
}
 
BGP_DMASend_Phase2::~BGP_DMASend_Phase2() {
    if (target_hosts_phase2_) {
        delete[] target_hosts_phase2_;
    }
}
 
// helps debugging when core dump since otherwise cannot tell where
// BGP_DMASend::send fails
#if 0
static void mymulticast(DCMF_Multicast_t* arg) {
    DCMF_Multicast(arg);
}
static void myrestart(DCMF_Request_t* arg) {
    DCMF_Restart(arg);
}
#endif
 
void BGP_DMASend::send(int gid, double t) {
    unsigned long long tb = DCMFTIMEBASE;
    if (NTARGET_HOSTS_PHASE1) {
        spk_.gid = gid;
        spk_.spiketime = t;
#if BGP_INTERVAL == 2
        bgp_receive_buffer[next_rbuf]->nsend_ += ntarget_hosts_;
        bgp_receive_buffer[next_rbuf]->nsend_cell_ += 1;
        if (next_rbuf == 1) {
            spk_.gid = ~spk_.gid;
        }
#else
        bgp_receive_buffer[0]->nsend_ += ntarget_hosts_;
        bgp_receive_buffer[0]->nsend_cell_ += 1;
#endif
        nsend_ += 1;
#if BGPDMA & 1
        if (use_bgpdma_) {
            nrnmpi_bgp_multisend(&spk_, NTARGET_HOSTS_PHASE1, target_hosts_);
        }
#endif
    }
    dmasend_time_ += DCMFTIMEBASE - tb;
}
 
void BGP_DMASend_Phase2::send_phase2(int gid, double t, BGP_ReceiveBuffer* rb) {
    unsigned long long tb = DCMFTIMEBASE;
    if (ntarget_hosts_phase2_) {
        spk_.gid = gid;
        spk_.spiketime = t;
#if BGP_INTERVAL == 2
        if (rb->index_ == 1) {
            spk_.gid = ~spk_.gid;
        }
#endif
        rb->phase2_nsend_cell_ += 1;
        rb->phase2_nsend_ += ntarget_hosts_phase2_;
#if BGPDMA & 1
        if (use_bgpdma_) {
            nrnmpi_bgp_multisend(&spk_, ntarget_hosts_phase2_, target_hosts_phase2_);
        }
#endif
    }
    dmasend_time_ += DCMFTIMEBASE - tb;
}
 
void bgp_dma_receive(NrnThread* nt) {
    //  nrn_spike_exchange();
    assert(nt == nrn_threads);
    TBUF double w1, w2;
    int ncons = 0;
    int& s = bgp_receive_buffer[current_rbuf]->nsend_;
    int& r = bgp_receive_buffer[current_rbuf]->nrecv_;
#if ENQUEUE == 2 && TBUFSIZE
    unsigned long tfind, tsend;
#endif
    w1 = nrnmpi_wtime();
#if BGPDMA & 1
    if (use_bgpdma_) {
        nrnbgp_messager_advance();
        TBUF
#if ENQUEUE == 2 && TBUFSIZE
            // want the overlap with computation, not conserve
            tfind = enq2_find_time_;
        tsend = enq2_enqueue_time_ - enq2_find_time_;
#endif
#if TBUFSIZE
        nrnmpi_barrier();
#endif
        TBUF while (nrnmpi_bgp_conserve(s, r) != 0) {
            nrnbgp_messager_advance();
            ++ncons;
        }
        TBUF
    }
#endif
    w1 = nrnmpi_wtime() - w1;
    w2 = nrnmpi_wtime();
#if TBUFSIZE
    tbuf_[itbuf_++] = (unsigned long) ncons;
    tbuf_[itbuf_++] = (unsigned long) bgp_receive_buffer[current_rbuf]->nsend_cell_;
    tbuf_[itbuf_++] = (unsigned long) s;
    tbuf_[itbuf_++] = (unsigned long) r;
    tbuf_[itbuf_++] = (unsigned long) dmasend_time_;
    if (use_phase2_) {
        tbuf_[itbuf_++] = (unsigned long) bgp_receive_buffer[current_rbuf]->phase2_nsend_cell_;
        tbuf_[itbuf_++] = (unsigned long) bgp_receive_buffer[current_rbuf]->phase2_nsend_;
    }
#endif
#if (BGPMDA & 2) && MAXNCONS
    if (ncons > MAXNCONS) {
        ncons = MAXNCONS;
    }
    ++xtra_cons_hist_[ncons];
#endif  // MAXNCONS
#if ENQUEUE == 0
    bgp_receive_buffer[current_rbuf]->enqueue();
#endif
#if ENQUEUE == 1
    bgp_receive_buffer[current_rbuf]->enqueue1();
    TBUF bgp_receive_buffer[current_rbuf]->enqueue2();
#endif
#if ENQUEUE == 2
    bgp_receive_buffer[current_rbuf]->enqueue();
    s = r = bgp_receive_buffer[current_rbuf]->nsend_cell_ = 0;
    bgp_receive_buffer[current_rbuf]->phase2_nsend_cell_ = 0;
    bgp_receive_buffer[current_rbuf]->phase2_nsend_ = 0;
    enq2_find_time_ = 0;
    enq2_enqueue_time_ = 0;
#if TBUFSIZE
    tbuf_[itbuf_++] = tfind;
    tbuf_[itbuf_++] = tsend;
#endif
#endif  // ENQUEUE == 2
    wt1_ = nrnmpi_wtime() - w2;
    wt_ = w1;
#if BGP_INTERVAL == 2
    // printf("%d reverse buffers %g\n", nrnmpi_myid, t);
    if (n_bgp_interval == 2) {
        current_rbuf = next_rbuf;
        next_rbuf = ((next_rbuf + 1) & 1);
    }
#endif
    TBUF
}
 
void bgp_dma_send(PreSyn* ps, double t) {
#if 0
    if (nrn_use_localgid_) {
        nrn_outputevent(ps->localgid_, t);
    }else{
        nrn2ncs_outputevent(ps->output_index_, t);
    }
#endif
    if (ps->bgp.dma_send_)
        ps->bgp.dma_send_->send(ps->output_index_, t);
}
 
void bgpdma_cleanup_presyn(PreSyn* ps) {
    if (ps && ps->bgp.dma_send_) {
        if (ps->output_index_ >= 0) {
            delete ps->bgp.dma_send_;
            ps->bgp.dma_send_ = 0;
        }
        if (ps->output_index_ < 0) {
            delete ps->bgp.dma_send_phase2_;
            ps->bgp.dma_send_phase2_ = 0;
        }
    }
}
 
static void bgpdma_cleanup() {
    nrntimeout_call = 0;
    for (const auto& iter: gid2out_) {
        bgpdma_cleanup_presyn(iter.second);
    }
    for (const auto& iter: gid2in_) {
        bgpdma_cleanup_presyn(iter.second);
    }
    if (!use_bgpdma_ && bgp_receive_buffer[1]) {
        delete bgp_receive_buffer[0];
        bgp_receive_buffer[0] = NULL;
    }
#if BGP_INTERVAL == 2
    if ((!use_bgpdma_ || n_bgp_interval != 2) && bgp_receive_buffer[1]) {
        delete bgp_receive_buffer[1];
        bgp_receive_buffer[1] = NULL;
    }
#endif
}
 
#ifndef BGPTIMEOUT
#define BGPTIMEOUT 0
#endif
 
#if BGPTIMEOUT
static void bgptimeout() {
    printf("%d timeout %d %d %d\n",
           nrnmpi_myid,
           current_rbuf,
           bgp_receive_buffer[current_rbuf]->nsend_,
           bgp_receive_buffer[current_rbuf]->nrecv_);
}
#endif
 
#if WORK_AROUND_RECORD_BUG
static void ensure_ntarget_gt_3(BGP_DMASend* bs) {
    // work around for bug in RecordReplay
    if (bs->ntarget_hosts_ > 3) {
        return;
    }
    int nold = bs->ntarget_hosts_;
    int* old = bs->target_hosts_;
    bs->target_hosts_ = new int[4];
    for (int i = 0; i < nold; ++i) {
        bs->target_hosts_[i] = old[i];
    }
    delete[] old;
    int h = (nrnmpi_myid + 4) % nrnmpi_numprocs;
    while (bs->ntarget_hosts_ < 4) {
        int b = 0;
        for (int i = 0; i < bs->ntarget_hosts_; ++i) {
            if (h == bs->target_hosts_[i]) {
                h = (h + 4) % nrnmpi_numprocs;
                b = 1;
                break;
            }
        }
        if (b == 0) {
            bs->target_hosts_[bs->ntarget_hosts_++] = h;
            h = (h + 1) % nrnmpi_numprocs;
        }
    }
    bs->ntarget_hosts_phase1_ = bs->ntarget_hosts_;
}
#endif
 
#define FASTSETUP 1
#if FASTSETUP
#include "bgpdmasetup.cpp"
#endif
 
void bgp_dma_setup() {
    bgpdma_cleanup();
    if (!use_bgpdma_) {
        return;
    }
    // not sure this is useful for debugging when stuck in a collective.
#if BGPTIMEOUT
    nrntimeout_call = bgptimeout;
#endif
    nrnmpi_bgp_comm();
    // if (nrnmpi_myid == 0) printf("bgp_dma_setup()\n");
    // although we only care about the set of hosts that gid2out_
    // sends spikes to (source centric). We do not want to send
    // the entire list of gid2in (which may be 10000 times larger
    // than gid2out) from every machine to every machine.
    // so we accomplish the task in two phases the first of which
    // involves allgather with a total receive buffer size of number
    // of cells (even that is too large and we will split it up
    // into chunks). And the second, an
    // allreduce with receive buffer size of number of hosts.
    max_ntarget_host = 0;
    max_multisend_targets = 0;
 
#if FASTSETUP
    // completely new algorithm does one and two phase.
    setup_presyn_dma_lists();
#else  // obsolete
    // see 672:544c61a730ec
#error "FASTSETUP required"
#endif  // obsolete (not FASTSETUP)
 
    if (!bgp_receive_buffer[0]) {
        bgp_receive_buffer[0] = new BGP_ReceiveBuffer();
    }
#if BGP_INTERVAL == 2
    if (n_bgp_interval == 2 && !bgp_receive_buffer[1]) {
        bgp_receive_buffer[1] = new BGP_ReceiveBuffer();
    }
#endif
}
 
#ifdef USENCS
 
// give me data on which gids of this node send out APs
int ncs_bgp_sending_info(int** sendlist2build) {
    int nsrcgid = 0;
    for (const auto& iter: gid2out_) {
        if (iter.second->output_index_ >= 0) {
            ++nsrcgid;
        }
    }
 
    *sendlist2build = nsrcgid ? new int[nsrcgid] : 0;
    int i = 0;
    for (const auto& iter: gid2out) {
        PreSyn* ps = iter.second;
        if (ps->output_index_ >= 0) {
            (*sendlist2build)[i] = ps->gid_;
            ++i;
        }
    }
}
}
 
// printf( "Node %d sees first hand %d gids send\n", nrnmpi_myid, nsrcgid );
return nsrcgid;
}
 
 
// function to access the sending information of a presyn
int ncs_bgp_target_hosts(int gid, int** targetnodes) {
    auto iter = gid2out_->find(gid);
    nrn_assert(iter != gid2out_.end());
    PreSyn* ps = iter->second;
    if (ps->bgp.dma_send_) {
        (*targetnodes) = ps->bgp.dma_send_->ntarget_hosts_
                             ? new int[ps->bgp.dma_send_->ntarget_hosts_]
                             : 0;
        return ps->bgp.dma_send_->ntarget_hosts_;
    }
 
    return 0;
}
 
// iterate over gid2in_ just so I can see what is in there
int ncs_bgp_target_info(int** presyngids) {
    //(*presyngids) = 0;
 
    int i, nsrcgid;
    nsrcgid = 0;
 
    // some target PreSyns may not have any input
    // so initialize all to -1
    for (const auto& iter: gid2in_) {
        PreSyn* ps = iter.second;
        assert(ps->output_index_ < 0);
        if (ps->bgp.srchost_ != -1) {  // has input
            ++nsrcgid;
            // printf( "Node %d: Presyn for gid %d has src %d\n", nrnmpi_myid, ps->gid_,
            // ps->bgp.srchost_ );
        }
    }
 
    (*presyngids) = nsrcgid ? new int[nsrcgid] : 0;
 
    i = 0;
    for (const auto& iter: gid2in_) {
        PreSyn* ps = iter.second;
        assert(ps->output_index_ < 0);
        if (ps->bgp.srchost_ != -1) {  // has input
            (*presyngids)[i] = ps->gid_;
            ++i;
        }
    }
 
    return nsrcgid;
}
 
int ncs_bgp_mindelays(int** srchost, double** delays) {
    int i, nsrcgid = 0;
 
    for (const auto& iter: gid2in_) {
        PreSyn* ps = iter.second;
        assert(ps->output_index_ < 0);
        if (ps->bgp.srchost_ != -1) {  // has input
            ++nsrcgid;
        }
    }
 
    (*delays) = nsrcgid ? new double[nsrcgid] : 0;
    (*srchost) = nsrcgid ? new int[nsrcgid] : 0;
 
    i = 0;
    for (const auto& iter: gid2in_) {
        PreSyn* ps = iter.second;
        assert(ps->output_index_ < 0);
        if (ps->bgp.srchost_ != -1) {  // has input
            (*delays)[i] = ps->mindelay();
            (*srchost)[i] = ps->bgp.srchost_;
            ++i;
        }
    }
 
    return nsrcgid;
}
 
#endif  // USENCS