Reduce

int MPI_Reduce(void *sbuf, void *rbuf, int count, MPI_Datatype dtype, 
               MPI_op op, int root, MPI_Comm comm)
{
  if (op is commutative) {
    return(reduce_commutative(sbuf, rbuf, count, dtype, op, root, comm));
  } else {
    return(reduce_noncommutative(sbuf, rbuf, count, dtype, op, root, comm));
  }
}

int reduce_commutative(void *sbuf, void *rbuf, int count, MPI_Datatype dtype, 
                       MPI_op op, int root, MPI_Comm comm)
{
  MPI_Status status;
  int myrank, mynum, rootnum, vnum, dim, hibit, mask;
  int peer, nmasters;
  void *tmpbuf, *redbuf;

  nmasters = num_masters(comm);
  MPI_Comm_rank(comm, &myrank);

  /* local phase */
  if (are_local(myrank, root, comm)) {
    perform a local reduction to root's rbuf;
  } else {
    perform a local reduction to a temporary buffer redbuf in
      the local master;
  }

  /* global phase */
  if ((myrank == root) || (is_master(myrank, comm) 
      && !are_local(myrank, root, comm))) {
  
    allocate a temporary buffer tmpbuf large enough for
      count copies of dtype;

    if (nmasters <= MAXLINEARREDUCE) {
      /* linear reduction to root */
      if (rank == root) {
        for (i = 0; i < nmasters; i++) {
          if (i == local_master_num(root, comm)) continue;
          MPI_Recv(tmpbuf, count, dtype, master_rank(i, comm), 
                   IMPI_REDUCE_TAG, comm, &status);
        }

        call reduction function op on tmpbuf (invec) 
          and rbuf (inoutvec);
      } else {
        MPI_Send(redbuf, count, dtype, root, IMPI_REDUCE_TAG, comm); 
      }
    } else {
      /* tree reduction to root */
      mynum = local_master_num(myrank, comm);
      rootnum = master_num(root, comm);
      vnum = (mynum - rootnum + nmasters) % nmasters;
      dim = cubedim(nmasters);
      hibit = hibit(vnum, dim);
      --dim;

      /* loop over cube dimensions */
      for (i = 0, mask = 1; i < dim; ++i, mask <<= 1) {
        /* a high-proc sends to low-proc and stops */
        if (vnum & mask) {
          peer = master_rank(((vnum & ~mask) + rootnum) % nmasters, comm);
          if (are_local(peer, root, comm)) {
            peer = root;
          }

          MPI_Send(redbuf, count, dtype, peer, IMPI_REDUCE_TAG, comm);
          break;
        }

        /* a low-proc receives, reduces, and moves 
         * to a higher dimension */
        else {
          peer = vnum | mask;
          if (peer >= nmasters) continue;
          peer = master_rank((peer + rootnum) % nmasters, comm);

          MPI_Recv(tmpbuf, count, dtype, peer, 
                   IMPI_REDUCE_TAG, comm, &status);

          if (myrank == root) {
            call reduction function op on tmpbuf (invec) 
              and redbuf (inoutvec);
          } else {
            call reduction function op on tmpbuf (invec) 
              and rbuf (inoutvec);
          }
        }
      }
    }
  }

  free all temporary buffers;
  return(MPI_SUCCESS);
}

int reduce_noncommutative(void *sbuf, void *rbuf, int count, 
                MPI_Datatype dtype, MPI_op op, int root, MPI_Comm comm)
{
  MPI_Status status;
  int i, myrank, nprocs;
  void *inbuf, *tmpbuf;

  MPI_Comm_size(comm, &nprocs);
  MPI_Comm_root(comm, &myrank);

  if (myrank != root) {
    return(MPI_Send(sbuf, count, dtype, root, IMPI_REDUCE_TAG, comm);
  }

  if (nprocs > 1) {
    create a temporary buffer tmpbuf large enough for count dtypes;
  }

  if (myrank == (nprocs - 1)) {
    MPI_Sendrecv(sbuf, count, dtype, myrank, IMPI_REDUCE_TAG,
                 rbuf, count, dtype, myrank, IMPI_REDUCE_TAG, comm, &status);
  } else {
    MPI_Recv(rbuf, count, dtype, nprocs - 1, IMPI_REDUCE_TAG, comm, &status);
  }

  for (i = nprocs - 2; i >= 0; --i) {
    if (myrank == i) {
       inbuf = sbuf;
    } else {
      MPI_Recv(tmpbuf, count, dtype, i, IMPI_REDUCE_TAG, comm, &status);
      inbuf = tmpbuf;
    }

    call reduction function op on inbuf (invec) and rbuf (inoutvec);
  }

  free all temporary buffers;
  return(MPI_SUCCESS);
}